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Abstract 


The objective of this project has been to improve the modeling of interactions between large 
Antarctic ice shelves and adjacent regions of the Southern Ocean. Our larger goal is to gain a 
better understanding of the extent to which the ocean controls ice shelf attrition, thereby 
influencing the size and dynamics of the Antarctic Ice Sheet. Melting and freezing under ice 
shelves also impacts seawater properties, regional upwelling and sinking and the larger-scale 
ocean circulation. Modifying an isopycnal coordinate general circulation model for use in 
sub-ice shelf cavities, we found that the abrupt change in water column thickness at an ice 
shelf front does not form a strong barrier to buoyancy-driven circulation across the front. 
Outflow along the ice shelf base, driven by melting of the thickest ice, is balanced by deep 
inflow. Substantial effort was focused on the Filch ner-Ronne cavity, where other models have 
been applied and time-series records are available from instruments suspended beneath the 
ice. A model comparison indicated that observed changes in the production of High Salinity 
Shelf Water could have a major impact on circulation within the cavity. This water propagates 
into the cavity with an asymmetric seasonal signal that has similar phasing and shape in the 
model and observations, and can be related to winter production at the sea surface. Even 
remote parts of the sub-ice shelf cavity are impacted by external forcing on sub-annual time 
scales. This shows that cavity circulations and products, and therefore cavity shape, will 
respond to interannual variability in sea ice production and longer-term climate change. The 
isopycnal model gives generally lower net melt rates than have been obtained from other 
models and oceanographic data, perhaps due to its boundary layer formulation, or the lack of 
tidal forcing. Work continues on a manuscript describing the Ross cavity results. 


Project Description 

We modified the Miami Isopycnal Coordinate Ocean Model (MICOM) for a domain that 
included the waters beneath an ice shelf and those occupying a portion of the open ocean. A 
floating ice shelf was incorporated as an upper boundary in the southern part of the domain. 
To initialize the model ocean, an extrapolation scheme filled the sub-ice cavity with waters 
similar to those observed during summer expeditions at the ice front. Simulations showed that 
initial properties were flushed out within a few years. Since year-round sea surface 
temperature and salinity observations do not exist over the model domain, we restored to 
summer measurements and to the sea surface freezing point in the open ocean during winter, 
adding sufficient salt during this period to generate the observed High Salinity Shelf Water. 

The sub-ice-shelf part of the domain is isolated from the atmosphere, but influenced by the 
pressure of the floating ice, and by phase changes that occur at the ice shelf base. Zero surface 
buoyancy forcing indicated that the introduction of a static pressure induced negligible motion 
in the underlying water. Nonzero surface buoyancy forcing produces a cyclonic circulation 
beneath the ice shelf. Following initial idealized simulations, realistic geometries were 
implemented for the Ross and Filchner-Ronne Ice Shelves and their underlying cavities. The 
numerical experiments included sensitivity studies to assess the fidelity of the simulations, 
and to check for natural internal oscillations and variability arising from external forcing. 
Model tracers were used to gain insight into system behavior. 



A key feature of the model is that vertical layering is based on density, with depth becoming a 
prognostic variable. This approach approximates the behavior of the real ocean, where 
advection and diffusion occur more readily along than across neutral or isopycnal surfaces. 
The model also has to accommodate strong buoyancy forcing, requiring that density be 
allowed to evolve freely in the upper layer. This variable density mixed layer applies both to 
the open ocean and ice cavity, requiring proper resolving of the transition between the two 
domains. Models of ocean circulation beneath ice shelves are driven in part by the heat and 
freshwater fluxes associated with phase changes at the ice-ocean boundary. One of the main 
differences between models is the treatment of turbulent transfer within the oceanic boundary 
layer. A comparison of alternative formulations gave results within -30% of observations 
beneath sea ice. However, the use of bulk transfer coefficients that specified turbulence can 
give melt/freeze rates that differ by as much as 5x. With our chosen formulation, net melting 
tends to be one half to one third of prior estimates. We have yet to include a tidal component, 
which can represent up to 20 cm/s mean flow in other models. 


Simulations of the Ross Ice Shelf cavity show that inflow is strongest on the western side, as 
would be expected from the denser shelf w'ater in that sector, with a relatively strong outflow 
near -175W. Several interior gyres occur in the vicinity of ice and sea floor topographic 
features, similar to those that have been reported beneath the other large ice shelves. The age 
of water at 450m is slightly less that the mean age estimated for the cavity from tracer studies 
in the open ocean. On a typical simulation, basal melting under the Ross Ice Shelf ranged 
from -0.3 to +0.3 m/yr, with melt patches resulting from peaks in the ice thickness, which are 
not altered during the model run. Such features might disappear with a higher-resolution ice 
thickness grid, but they also suggest that net melting is likely to be dominated by thicker ice 
where it moves off the grounding line. The annual cycle of monthly averaged net melting for 
the Ross cavity shows a peak in the late winter/early spring, confirming a seasonal cycle that 
appeared earlier in a 2-D channel model around Roosevelt Island. 

As the Antarctic coastline also contains smaller ice shelves and floating glaciers of various 
sizes and different hydrographic environments, we carried out a simulation for the Pine Island 
Glacier (PIG) which is located in Pine Island Bay in the Amundsen Sea. Initial estimates of 
the shape and size of the PIG cavity, coupled with ocean measurements taken during a 1994 
field experiment, resulted in a modeled net melt rate of -14 m/yr. This is -2 orders of 
magnitude higher than simulated for the Ross cavity with the same model, but is consistent 
with independent estimates that range from 10-24 m/yr for the PIG cavity. Sensitivity studies 
confirmed earlier inferences that the high melt rate was the combined result of the ‘warm’ 

(>+ 1 C) Circumpolar Deep Water that penetrates the cavity, and the relatively steep slope of 
the ice shelf base. A general review of the Pine Island Glacier system did not reveal clear 
evidence for retreat or collapse over the last few decades. 


Products 

This report incorporates the primary publications that have resulted in full or in part from 
NASA support of this project. Those papers and related manuscripts submitted or in 
preparation are listed in the accompanying bibliography. Presentations of the ongoing work, 
listed following the bibliography, were made at an AGU meeting in San Francisco, a West 
Antarctic Ice Sheet Workshop in Sterling VA, a MICOM Workshop in Miami, FL, a 
Chapman Conference in Orono, ME, a Gordon Research Conference in Ventura, CA, a Ross 
Sea Oceanography Conference in Ischia, Italy, and at the European Geophysical Society in 
Nice, France. PI Holland also traveled to other sites for talks and discussions, and co-PI 
Jenkins was supported for travel-related expenses during time spent at Lamont working on 
this project. 



About midway through the grant, PI Holland assumed a faculty position at the Courant 
Institute of Mathematics of New York University. This necessitated some realignment of 
expenses, and travel between NYU and Lamont, where Holland was granted an adjunct 
appointment to facilitate continued work on the project. 
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Introduction 

The cavities beneath the floating ice shelves of Antarctica are regions where the waters of the 
Southern Ocean come into direct contact with the ice sheet. Exchanges of mass and energy at 
the ice/ocean interface effect a transformation of water masses flowing into the cavity and 
contribute to the wastage side of Antarctica’s mass balance equation. How sensitive the 
conditions in the sub-ice cavities are to external changes is a question of some importance, as the 
past and future evolution of both the ice sheet and the surrounding seas are at least partially linked 
to the answer. Our understanding of the processes that drive sub-ice-shelf circulation owes muc 
to the application of computer models of varying sophistication to the problem (Williams et a 
in press) Recently Grosfeld et al. ( 1 997) presented results from the first model to include both 
an open ocean and a sub-ice cavity within its domain. They concluded that the wind-forced 
circulation north of the ice shelf and the buoyancy-forced circulaUon beneath the ice shelf were 
largely separate, with the degree of communication being topographically controlled. Here we 
describe some of the preliminary results we have obtained applying an isopycnic coordinate ocean 

model to a similar domain. 


MICOM 

We have used the Miami Isopycnic Coordinate Ocean Model (MICOM; Bleck et al, 1992) tor 
this investigation. The key feature of this model is that the vertical discretisation of the ; domain 
is based on density instead of the more usual depth. In the isopycnic formulation depth becomes 
a prognostic variable. Whereas conventional ocean models calculate the water density at a fixed 
number of discrete depths, MICOM computes the depth at which a fixed number of discrete 
densities occur. In practice, the integration is performed layer-wise; that is, the ocean is dm e 
into a number of layers of constant-density water, and the model computes the thickness of he 
layers at each point of the horizontal grid. The advantage of this approach is that it mimics the 
behaviour of the real ocean very well. Advection and diffusion in the ocean occur readily along 
neutral surfaces (which are approximately isopycnal surfaces), but much less easily across them. 
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The disadvantage is that if a given water density does not occur at certain points in the domain, 
that particular layer will have regions of zero thickness. A sophisticated numerical scheme is 
required to handle this situation without generating negative layer thicknesses. 

If an isopycnic model is to accept buoyancy forcing, the upper layer must have a freely-evolving 
density. A mixed layer formulation based on the work of Gaspar (1988) is used for the surface 
layer in MICOM. The entrainment/detrainmeni algorithm allows water to be cycled from one 
isopycnic layer, through the mixed layer, where its properties are modified by the surface 
boundary conditions, into another isopycnic layer. This procedure mimics the manner in which 
water mass transformations occur in nature, but herein lies the greatest hurdle to be overcome in 
the application of MICOM to a sub-ice- shelf cavity. The boundary layer adjacent to the ice shelf 
base is also a site of buoyancy forcing and water mass modification, so the mixed layer must 
extend under the ice shelf. We have made extensive changes to the model code to allow this to 
happen. 


Testing the new code 

We have run a number of simulations to test the stability and accuracy of our new code: 

(i) We have set up an ocean with horizontal isopycnal surfaces in front of and beneath an ice 
shelf to check that the water can remain stationary if not forced. 

(ii) We have placed an isolated patch of light fluid beneath the ice shelf to check that it can 
flow up the ice shelf base and escape from the cavity. 

(iii) We have applied boundary conditions at the ice shelf base that freshen and cool the mixed 
layer to check that its response is to flow up the ice shelf base and exit the cavity. 

All of these tests were applied successfully, both with and without Earth rotation. 


Preliminary results 

Our first results with the full model have teen obtained on a 10° latitude by 15° longitude domain, 
with a horizontal resolution of 20-35 km and 10 layers in the vertical. The seabed was flat and 
1000 m deep everywhere, while an ice shelf covered the southern half of the domain. The ice 
thickness was zonally invariant, but increased from 0 m at the ice front to 300 m at the southern 
edge of the domain. The water column was initially 0°C everywhere, with salinity ranging from 
34.4-34.8 psu, and the model was spun up for five years. 

The basal freezing rate (negative numbers indicate melting) at the base of the ice shelf is illustrated 
in Figure 1 . No freezing occurs because the inflowing water is almost 2°C above the surface 
freezing point and peak melt rates, in excess of 3.5 m yr' 1 , occur in the south-eastern corner of 
the domain. The stream function for the vertically integrated transport shows a single, western- 
intensified gyre occupying most of the cavity (Figure 2). It extends beyond the ice front, which 
crosses the domain at 75°S, but the currents to the north of the ice shelf are considerably weaker. 
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Figure 1: Basal freezing rate (m yr' 1 ) at the ice shelf base. 


This circulation pattern is broadly similar to that found by Determann and Gerdes (1994) and 
Grosfeld et al. (1997). However, we need to look a little closer at the model results in order to 
be able to explain the pattern of basal melting. If the entire water column were following this gyre 
we would expect strongest melting immediately to the east of the gyre centre, where the inflowing 
currents are strongest. The most active layer beneath the ice shelf is the mixed layer, which reacts 
to the combination of buoyancy and Coriolis forcing by flowing from east to west across the ice 
shelf base, with a small component to the north generated by the frictional drag (Figure 3). This 
leads to strong divergence along the eastern boundary and weak divergence along the southern 
boundary. The resulting upwelling of warm water and entrainment into the mixed layer gives rise 
to the pattern of melting and sustains the motion. The mixed layer behaves in the manner we 
would expect of a buoyant, two dimensional plume under the influence of rotation. This suggests 
that the differences between the circulation produced by one and two dimensional models (eg. 
Jenkins, 1991; Hellmer and Olbers, 1989) may not be as great as suggested by Williams et al. (in 
press). 





Figure 2: Stream function for the vertically integrated transport (Sv). 


More results from this model run can be found on the web site: 

http://figgy.ldeo.columbia.edu/~holland/html_float_sheliyidealized/cxpt_002/expt_descriptor.html 


Future work 

These preliminary results suggest that our modified version of M1C0M can be a useiul tool in 
understanding the intricacies of sub-ice-shelf circulation. In future we hope to apply the model 
to domains more reminiscent of the southern Weddell Sea, and eventually to use real topography. 
As M1C0M is a free surface model we also hope to be able to include the effects of tidal forcing 

at a later date. 
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ABSTRACT 

Models of ocean circulation beneath ice shelves are driven primarily by the heat and freshwater fluxes that 
are associated with phase changes at the ice-ocean boundary. Their behavior is therefore closely linked to the 
mathematical description of the interaction between ice and ocean that is included in the code. An hierarchy of 
formulations that could be used to describe this interaction is presented. The main difference between them is 
the treatment of turbulent transfer within the oceanic boundary layer. The computed response to various levels 
of thermal driving and turbulent agitation in the mixed layer is discussed, as is the effect of various treatments 
of the conductive heat flux into the ice shelf. The performance of the different formulations that have been used 
in models of sub-ice-shelf circulation is assessed in comparison with observations of the turbulent heat flux 
beneath sea ice. Formulations that include an explicit parameterization of the oceanic boundary layer give results 
that lie within about 30% of observation. Formulations that use constant bulk transfer coefficients entail a definite 
assumption about the level of turbulence in the water column and give melt/freeze rates that vary by a factor 
of 5, implying very different forcing on the respective ocean models. 


1. Introduction 

The continental shelf seas surrounding Antarctica 
most frequently attract attention because they are the 
source regions of Antarctic Bottom Water. It is com- 
monly assumed that atmospheric forcing of the ocean 
and ice cover is the primary driving mechanism behind 
the deep convection that occurs over the continental 
slope (e.g., Gill 1973). However, poleward of the shelf 
break, 40% of the sea surface is covered by floating ice 
shelves, which range in thickness from 100 to 2000 m 
and therefore completely isolate the ocean from the at- 
mosphere. Circulation beneath the ice shelves and the 
associated meltwater input have a profound impact on 
shelf water properties (Foldvik et al. 1985; Jacobs et al. 
1985; Fahrbach et al. 1994; Hellmer et al. 1999). Togg- 
weiler and Samuels (1995) suggest that up to 75% of 
all the ocean's deep waters may retain a signature of 
this meltwater input. 

The interaction between ice shelves and the ocean is 
thus a potentially important element of the climate sys- 
tem, and in recent years numerical models have been 
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used to evaluate the key processes operating in sub-ice- 
shelf cavities (Williams et al. 1999). Upper boundary 
conditions derived from a thermodynamic model of the 
ice-ocean interaction have been applied to ocean models 
of varying sophistication. Dynamic models of the ice 
shelf itself have not been included, so the ice-ocean 
interface has been treated as a fixed boundary. The dis- 
parity of timescales between the slowly flowing ice shelf 
and the relatively fast flowing waters beneath provides 
some justification for this approach. Although the spec- 
ification of the upper boundary conditions represents a 
computationally small and simple component of most 
sub-ice-shelf circulation models, it is of crucial impor- 
tance. A correct estimate of the surfaces fluxes is es- 
sential to a realistic simulation of the sub-ice-shelf cir- 
culation and to the utility of the results for estimating 
the mass balance of the ice shelves. 

In this paper we focus on the mathematical descrip- 
tion of the ice-ocean interaction. We present an hier- 
archy of models describing the heat and freshwater ex- 
change at and near the ice-ocean interface. Our aim is 
to compare the behavior of the differing upper boundary 
formulations that have been used in models of sub-ice- 
shelf circulation to date and to introduce a new for- 
mulation that closely follows the work of McPhee et al. 
(1987). For a related comparative study, but in the con- 
text of sea ice-ocean coupling, the reader is referred to 


© 1999 American Meteorological Society 



1788 


JOURNAL OF PHYSICAL OCEANOGRAPHY 


Volume 29 


work by Holland et al. (1997). We also compare the 
results of the various models with recent observations 
of heat flux in the turbulent boundary layer beneath sea 
ice. 

The fundamental assumption in all the models is that 
phase changes occur in thermodynamic equilibrium so 
that the temperature and salinity at the ice-ocean in- 
terface are always related by an expression for the freez- 
ing point at the appropriate depth. The problem becomes 
one of calculating the heat and freshwater fluxes that 
result from deviations in the far-field ocean properties 
from freezing point conditions. Our treatments of the 
processes occurring in the oceanic boundary layer differ 
only in the sophistication with which turbulent diffusion 
of heat and salt is modeled. On the ice side of the in- 
terface the flow is laminar and only the diffusion of heat 
need be considered because salt cannot diffuse through 
the solid ice matrix. The problem is that diagnosis of 
the temperature gradient at the ice shelf base requires 
solution of the equation for heat advection and diffusion 
throughout the ice shelf. As the ice flow is unknown, 
we require a reduced form of this equation that is trac- 
table but captures the essential features of the full so- 
lution. 

In common with most models of sub-ice-shelf cir- 
culation, we consider all phase changes to occur at the 
ice-ocean boundary, regardless of whether the far-field 
ocean conditions are above or below freezing. We do 
not consider the process of frazil ice growth in super- 
cooled parts of the water column, although observations 
and modeling (Oerter et al. 1992; Bombosch and Jenkins 
1995) suggest that deposition of suspended ice crystals 
is the dominant mechanism of basal growth beneath ice 
shelves. Thus, while it is instructive to intercompare the 
response of the different boundary formulations to su- 
percooling in the ocean, direct comparisons between 
models and observations are only valid for melting con- 
ditions. It is possible to treat the thermodynamics of 
frazil ice growth in a manner analogous to that outlined 
below for melting and freezing at a solid boundary, but 
the incorporation of frazil ice into an ocean model re- 
quires the addition of an ice conservation equation (Om- 
stedt and Svensson 1984; Mellor et al. 1986; Jenkins 
and Bombosch 1995). 

2. Thermodynamic models of ice-ocean interaction 

The objective of modeling the ice-ocean interaction 
is to obtain as realistic as possible a melt rate at the ice- 
shelf base. We now define and solve the necessary equa- 
tions to achieve this. The far-field quantities are the 
prescribed interior properties of the ice shelf and the 
properties of the upper layer or level of the ocean model. 
We are interested in determining the characteristics ex- 
actly at the ice-ocean interface where there are three 
physical constraints: the interface must be at the freezing 
point and both heat and salt must be conserved at the 
interface during any phase changes. This gives a system 


of up to three equations in up to three unknowns, name- 
ly, the interface temperature, salinity, and melt rate. 

To assist the discussion below, the relevant layers, 
temperatures, salinities, and heat and salt fluxes are 
shown schematically in Fig. 1. The ocean mixed layer 
has a temperature T M and salinity S M , which are not 
necessarily equal to the respective ice-ocean interface 
properties T B and S B . The gradients in temperature and 
salinity through the boundary layer drive heat and salt 
fluxes between the interface and the mixed layer. The 
temperature gradient in the ice at the base of the ice 
shelf drives a heat flux from the interface into the ice 
shelf, which has a surface temperature denoted by T s 
and a bulk salinity denoted by S f . 

a . Fundamental equations 

1) Freezing point dependence 

The freezing point of seawater is a weakly nonlinear 
function of salinity and a linear function of pressure 
(Millero 1978). This relationship between temperature 
and salinity at the ice-ocean interface will be one of 
three equations that will have to be solved simulta- 
neously, so it is simpler to work with a linearized ver- 
sion: 

T b = aS B + b + cp B , (1) 

where p B is the pressure at the interface. The values of 
the empirical constants a, b, and c are given in Table 1 
along with other constants and parameters used in this 
study. The formula is valid only in the salinity range 
4-40 psu and does not apply to pure freshwater. 


2) Heat conservation 

At the ice-ocean interface, the divergence of the heat 
flux balances the sink or source of latent heat caused 
by melting or freezing: 

Q] ~ Ql* = QL nr (2) 

The latent heat term i* given by 

Glatent — (3) 

where p M w B represents the mass of ice that is melted 
(w, > 0) or frozen (w g < 0) per unit time. The esti- 
mation of the diffusive heat fluxes is discussed in detail 
below. 

3) Salt conservation 

An equation analogous to (3) describes the salt flux 
required to maintain the boundary salinity at S B in the 
presence of the “freshwater” flux associated with melt- 
ing or freezing of ice having a salinity of S { : 

Qb nne ~ Pm W b($1 “ ^ B )* (4) 

This is balanced by the salt flux divergence at the in- 
terface: 
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Fig. 1. Schematic representation of (a) the heat and (b) the salt balance at the base of an ice 
shelf. The slope of the ice shelf base is greatly exaggerated for illustrative purposes. 


Q S , - Qlf = Q £«- ( 5 ) 

The diffusive flux of salt into the ice shelf, Q s , , is iden- 
tically zero and will be discussed no further, while the 
estimation of the diffusive flux through the oceanic 
boundary layer will be addressed below. Continental ice 
that melts from the base of an ice shelf has a salinity 
of zero, but when seawater freezes, brine is usually 
trapped within the forming ice, giving it a nonzero bulk 
salinity. Observations of marine ice found at the base 
of ice shelves with salinities of 0.025 psu (Eicken et al. 
1994) indicate that a very effective desalination process 
must operate. As further evidence, additional observa- 
tions (Oerter et al. 1992) show very low salinities of 


approximately 0.100 and suggest that to a good ap- 
proximation we can treat S, as zero always. 

b , Modeling the oceanic fluxes 

1) A ONE-EQUATION FORMULATION 

The simplest approach recognizes that whatever the 
details of heat and salt transfer through the oceanic 
boundary layer, the overall effect is to cause the upper 
layer of the ocean to relax toward the freezing point. If 
this relaxation is assumed to occur instantaneously, there 
is no distinction between interface and mixed layer 
properties, and the ice-ocean interaction can be de- 
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Table !. Model parameters and constants. 


Parameter 

Symbol 

Units 

Value 

Salinity coefficient of freezing equation 

a 

°C psu 1 

-5-73 X 10 2 

Constant coefficient of freezing equation 

b 

°C 

9.39 X 10 2 

Bernoulli numbers 

B n 

dimensionless 


Ocean surface buoyancy flux 

B„ 

m z s' 2 


Pressure coefficient of freezing equation 

c 

°C Pa 1 

-7.53 X 10 s 

Constants of integration 

G. c 2 

dimensionless 


Momentum exchange coefficient 

Cj 

dimensionless 

1.50 X 10 1 

Specific heat capacity ice shelf 

c pi 

Jkg 1 K -1 

2009.0 

Specific heat capacity mixed layer 

C P M 

J kg 1 K ! 

3974.0 

Coriolis parameter 

f 

s _l 

-1.00 X 10 4 

Gravitational acceleration 

S 

m $^ 2 

9.81 

Thickness of boundary layer 

h 

m 


Thickness of viscous sublayer 

h. 

m 

-0.001 

Thickness of ice shelf 

H, 

m 

-1000 

Thickness of mixed layer 

Hm 

m 

-10 

Von Kdrman’s constant 

k 

dimensionless 

0.40 

Latent heat fusion 

L, 

Jkg 1 

3.34 X 10 s 

Obukhov length 

L u 

m 


Solutions of characteristic polynomial 

m, m j, m 2 

dimensionless 


Nusselt number 

Nu 

dimensionless 

:»! 

Pressure at ice shelf base 

Pb 

Pa 

-1.0 X 10 7 

Prandtl number 

Pr 

dimensionless 

13.8 

Conductive heat flux through ice shelf 

QJ 

W m ? 


Latent heat at ice-ocean interface 

QLttnt 

W m- ? 


Diffusive heat flux in boundary layer 

Q T u 

W m 2 


Diffusive salt flux through ice shelf 

Q) 

psu m~‘ s _l 


Salt flux at ice-ocean interface 

Qlnn 

psu m _l s ’ 


Diffusive salt flux in boundary layer 

Qlt 

psu nr 1 s ’ 


Ratio of melt/freeze rates 

r 

dimensionless 


Critical flux Richardson number 

R, 

dimensionless 

0.20 

Schmidt number 

Sc 

dimensionless 

2432 

Bulk salinity of ice shelf 

S, 

psu 

0 

Salinity at ice-ocean interface 


psu 

Prognostic 

Salinity of mixed layer 

S u 

psu 

-34.5 

Time coordinate 

1 

s 


Temperature at ice shelf surface 

T s 

°C 

25.0 

Temperature at ice-ocean interface 

T, 

°C 

Prognostic 

Temperature of mixed layer 

T u 

°c 

1.85 

Thermal driving 

T, 

. °c 


Ice shelf flow velocity 

u, 

m s"' 


Ocean mixed layer velocity 


m s' 


Friction velocity ice-ocean 

U * 

m 


Melt rate at ice shelf base 


m s' 1 

Prognostic 

Vertical velocity of ice shelf 

*'/ 

ms 1 


Vertical geopotential coordinate 

z 

m 


Salinity contraction coefficient 

A, 

psu 1 


Thermal expansion coefficient 

Ar 

°c-‘ 


Molecular thermal conductivity ice shelf 

Kf 

m 2 s-' 

1.14 X 10 6 

Thermal diffusivity of mixed layer 
Molecular salt conductivity ice shelf 

Km 

K? 

m 2 

0.0 

Salt diffusivity mixed layer 

K« 

m 2 


Ice shelf reference density 

P, 

kg m 1 

920.0 

Ocean reference density 

Pm 

kg m 1 

1025 0 

Ratio of thermal driving to thermal forcing 

0 

dimensionless 


Thermal exchange velocity 

Jr 

m s~‘ 

-LOO X 10 4 

Salinity exchange velocity 

7s 

m s - ’ 

-5.05 X 10 7 

Kinematic viscosity of sea water 

V 

m 2 s' 1 

1.95 X 10 

Pec let number 

Y 

dimensionless 


Temperature gradient amplification factor 

n 

dimensionless 


Turbulent transfer parameter 

^Turfe 

dimensionless 


Thermal molecular transfer parameter 

VT 
1 Mtik 

dimensionless 


Salinity molecular transfer parameter 

ri* 

dimensionless 


Stability parameter 

% 

dimensionless 

<1 

Stability constant 

L 

dimensionless 

0.052 
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scribed completely using Eq. (1). Such a formulation 
has been widely used in large-scale ocean-atmosphere- 
sea ice models (Holland 1998) and was the basic as- 
sumption behind the earliest conceptual and numerical 
models of ice shelf-ocean interaction (Doake 1976; 
Robin 1979; MacAyeal 1985; Jenkins and Doake 1991). 
Despite its simplicity, the application of such a bound- 
ary condition to an ocean model may not be straight- 
forward. If the usual prognostic equations for temper- 
ature and salinity are solved everywhere, the derived 
values must subsequently be reset wherever the mixed 
layer is in contact with ice. The melt rate cannot be 
recovered from the boundary condition but is deter- 
mined from the change in temperature of the mixed 
layer. The derived rate is therefore a function of model 
time step and mixed layer thickness among other things. 
This is an undesirable feature if the aim is an accurate 
diagnosis of the melting and, since this formulation is 
not directly comparable to those presented below, we 
will discuss it no further. 

2) Two-equation formulations 

The next approach that we will discuss recognizes 
that the rate at which the mixed layer temperature re- 
laxes toward the freezing point is governed by the dif- 
fusion of heat through the oceanic boundary layer. Equa- 
tion (2) is introduced and estimates are made of the two 
heat fluxes that appear on the left-hand side. In their 
most general form they can be written 


the boundary layer means that the temperature profile 
is nonlinear and the diffusivity is variable, but we can 
parameterize these complications with the introduction 
of a Nusselt number, Nu, an empirical parameter having 
a value greater than 1 : 


/NukJ, 

Qm — Pm C pm\ h 


(T b - T„y 


( 9 ) 


dT, 

QJ = -P' c p' k *Yz 


and 


Qm — Pm c p m k m 


dz 


( 6 ) 


( 7 ) 


In the above equations k are thermal diffusivities ad- 
jacent to the ice-ocean interface, p are densities, and c p 
are specific heat capacities, while subscript M indicates 
mixed layer properties and / ice properties. In Eq. (6) 
the density, specific heat capacity, and thermal diffu- 
sivity may all be regarded as constant, and the problem 
becomes one of estimating the temperature gradient at 
the base of the ice shelf. To solve Eq. (7) we can treat 
the density and specific heat capacity as constant, but 
we need a suitable parameterization of the product of 
the diffusivity and temperature gradient near the ice 
shelf base. If the boundary layer were laminar, we would 
anticipate that the temperature would vary linearly be- 
tween the interface and mixed layer temperatures, in 
which case Eq. (7) could be written 


We will refer to the first quantity in brackets, which has 
dimensions of velocity, as a thermal exchange velocity 
y The simplest approach would be to choose a constant 
value for the exchange velocity but, recognizing that it 
is a result of turbulence in the mixed layer, a more 
realistic assumption is to make it a function of the fric- 
tion velocity. 

The two-equation formulation offers some advantag- 
es over the one-equation formulation. It is still very 
simple, but it includes a diagnosis of the melt rate. As- 
sociated heat and freshwater fluxes can then be applied 
to the ocean model in an identical manner to all other 
surface fluxes, and no special treatment of the mixed 
layer equations is required because of the presence of 
ice. However, it still lacks some realism, as to solve Eq. 
(1) it must be assumed that the interface salinity and 
the mixed layer salinity are identical. This implies in- 
finite salt diffusivity, whereas in reality we would expect 
salt to diffuse at the same rate as, or slower than, heat. 
Nevertheless, the only error implied by this assumption 
is the misdiagnosis of the interface temperature, and as 
this is a weak function of salinity, we might anticipate 
relatively small errors. McPhee (1992) and McPhee et 
al. (1999) show that this formulation produces heat flux- 
es that agree well with measurements made beneath sea 
ice having a wide range of roughness characteristics. 
An analogous formulation, but with a constant thermal 
exchange velocity and a constant, prescribed interface 
salinity, has been used in the sub-ice-shelf models of 
Determann and Gerdes (1994), Grosfeld et al. (1997), 
and Williams et al. (1998). 

3) Three-equation formulations 

The most sophisticated formulations make no prior 
assumptions about conditions at the interface and solve 
Eqs. (1), (2), and (5) using the known mixed layer and 
ice properties. By an analogous argument to that used 
in the derivation of Eq. (9) we can express the surviving 
term on the left-hand side of Eq. (5) as 


Q s m — ~PmJs(^b ^ m ). 


( 10 ) 


Qli = -P»c pM k t m 


T (Js - TJ 


( 8 ) 


where h is the boundary layer thickness. Turbulence in 


where y s is the salinity exchange velocity. Once again 
the exchange velocity can be either assumed constant 
or assigned a functional dependence on friction velocity. 
The former approach was followed by Hellmer and 01- 
bers (1989, 1991), Scheduikat and Olbers (1990), and 
Hellmer and Jacobs (1992). The latter approach has been 
developed in the sea-ice literature [see Gade (1993) for 
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a review] and adopted in the ice shelf-ocean interaction 
models of Jenkins (1991), Jenkins and Bombosch 
(1995), Hellmer and Jacobs (1995), and Hellmer et al, 
(1999). 


c. Parameterizing the transfer of heat and salt 

through the oceanic boundary layer 

The key to diagnosing a realistic melt rate from either 
the two- or the three-equation formulation lies in the 
choice of appropriate exchange velocities. In the case 
of the three-equation model the problem is complicated 
by the fact that the thermal and salinity diffusivities can 
only be assumed to be equal in the fully turbulent part 
of the boundary layer. Close to the ice-ocean interface, 
the eddy size and hence the turbulent diffusivity are 
suppressed. Where the suppression is great enough that 
molecular diffusion becomes the dominant transfer 
mechanism, heat will diffuse more rapidly than salt. As 
the exchange velocities need to account for all processes 
occurring within the boundary layer, y, will be smaller 
than y T . 

The role of molecular diffusion in governing the rate 
of heat and mass transfer within a thin, viscous sublayer 
adjacent to the ice-ocean boundary was recognized by 
Mellor et al. (1986). Subsequently, McPhee et al. (1987) 
and Steele et al. (1989) investigated boundary layer pa- 
rameterizations that explicitly included a viscous sub- 
layer. Jenkins (1991) used an analogous parameteriza- 
tion to calculate exchange velocities at the base of an 
ice shelf. Assuming the ice-ocean interface in this case 
to be hydraulically smooth leads to expressions of the 
form (Kader and Yaglom 1972) 

YT 2.12 In (u*hfv) + 12.5 Pr 2/3 - 9 
and 


7s 2.12 In (u*hlv) + 12.5 Sc 2 ' 3 - 9' 

The influence of the molecular sublayer is apparent in 
the inclusion of the molecular Prandtl number, Pr (the 
ratio of viscosity to thermal diffusivity), and the mo- 
lecular Schmidt number, Sc (the ratio of viscosity to 
salinity diffusivity), in the denominators. The kinematic 
viscosity of sea water is considered constant and is de- 
noted by the symbol v. The friction velocity w* is de- 
fined in terms of the shear stress at the ice-ocean in- 
terface, a simple parameterization of which involves a 
dimensionless drag coefficient c d and the velocity of the 
mixed layer U M , the ice being considered stationary: 

= cJU\r (13) 

A potentially important effect not accounted for in (1 1) 
and (12) is the impact of the buoyancy flux at the ice- 
ocean interface on turbulence within the boundary layer. 
A stabilizing buoyancy flux (i.e., melting) will suppress 


mixing, while a destabilizing buoyancy flux (i.e., freez- 
ing) will enhance mixing (McPhee 1994). We wish to 
investigate how the stability of the boundary layer might 
influence melt rates at the base of an ice shelf, so we 
follow McPhee et al. (1987) in expressing the transfer 
coefficients as 


7t.s = 


r -f vt.s 1 

1 Turb ' 1 Mote 


(14) 


where 


r = 1 . | 1 1 

Tu,b k { fh, I 2&1J, k 


(15) 


and 


TJL = 12.5(Pr, Sc) 273 - 6. (16) 

In Eq. (15) k is von Karmen’s constant, / is the Coriolis 
parameter, is a dimensionless constant, and h v is the 
thickness of the viscous sublayer. Values for the first 
three of these are given in Table 1, while we estimate 
the sublayer thickness to be (Tennekes and Lumley 
1972, p. 160) 


K = 5—. (17) 

In Eqs. ( 1 5)— ( 1 7) we have assumed the ice-ocean in- 
terface to be hydraulically smooth. The influence of the 
interfacial buoyancy flux is encapsulated in the stability 
parameter, introduced by McPhee (1981): 


* 7 * = 



fLoK) 


- 1/2 


(18) 


where R c is a critical flux Richardson number and L 0 is 
the Obukhov length. If the Obukhov length is negative 
(i.e., the buoyancy flux is destabilizing) the stability 
parameter is set to 1. We do not consider how desta- 
bilizing buoyancy fluxes influence freezing rates be- 
cause direct freezing to the ice shelf base is thought to 
be limited. The formation of frazil ice in the water col- 
umn, which we do not discuss here, will have a stabi- 
lizing effect on the boundary layer (Jenkins and Bom- 
bosch 1995). 


d . Modeling the heat flux into the ice shelf 


We now address the problem of estimating the basal 
temperature gradient in the ice shelf, required in Eq, 
(6). This requires solution of the heat transport equation 
in the ice shelf: 


— + U, • VT, = k[V 2 T,. (19) 

dt 

In practice, the solution of the full equation is not pos- 
sible unless the flow field within the ice shelf is known, 
so we consider reduced forms. 
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1) No ADVECTION, NO DIFFUSION 

The simplest of all approximations is that the ice shelf 
is a perfect insulator. With no diffusion into the ice shelf, 
the first term on the left-hand side of Eq. (2) is iden- 
tically zero. Such an approximation has been used by 
Determann and Gerdes (1994), Jenkins and Bombosch 
(1995), Grosfeld et al. (1997), and Williams et al. 
(1998), although it can be justified only if the conducted 
heat flux is always small compared to the latent heat 
term. 


2) No ADVECTION, VERTICAL DIFFUSION 
In this case (19) reduces to 

^ = 0 ( 20 ) 

dz 2 

for a steady state. The solution is a linear temperature 
profile throughout the thickness of the ice shelf, so the 
basal gradient can be expressed as 

dz 

where H l is the thickness of the ice shelf and T s is the 
surface temperature. Such an approximation has been 
used in the models of Hellmer and Olbers (1989, 1991), 
Scheduikat and Olbers (1990), and Hellmer and Jacobs 
(1992, 1995). This approach is common in sea-ice mod- 
els, but the thickness of ice shelves and the observed, 
highly nonlinear temperature profiles make it less sat- 
isfactory for modeling ice-shelf thermodynamics. 


(Ts - t b ) 
H, ’ 


(21) 


3) Constant vertical advection, vertical 

DIFFUSION 


The simplest way to allow a nonlinear temperature 
profile to develop is to allow for vertical advection with- 
in the ice shelf. We will assume that the vertical velocity 
is constant and equal to the basal melt/freeze rate and 
that the ice shelf is in a steady state. There is some 
justification for this approximation in that over the short 
timescales of interest here, that is, of years to decades, 
the ice shelves are believed to be in a relatively steady 
state. Clearly, such an approximation would be unre- 
alistic for consideration of sea ice thermodynamics be- 
cause of its thinness. This requires that all ice added or 
removed at the base is balanced by surface ablation or 
accumulation. In this case Eq. (19) reduces to the equa- 
tion used by Wexler (1960) (see also discussion by Pat- 
erson 1994, p. 204): 


where 


= 0 . 

dz 2 K] dz 


w, 



( 22 ) 


(23) 



Fig. 1 . Temperature-depth profiles through an ice shelf 1000 m 
thick calculated assuming a constant vertical velocity. Surface and 
basal’temperatures are -25° and -2°C, respectively. Vertical velocity 
(in m yr ') is given for each profile where positive labels indicate 

basal melting. 


We assume a solution of the form T(z) = which 
yields a quadratic characteristic polynomial of Eq. (22) 
having the two roots m, = -w,IkJ and m 2 = 0. The 
general solution of Eq. (22) is then 

T,(z ) = c t e m ' : + c 2 e mj; , (24) 

where c, and c 2 are constants to be determined by the 
boundary conditions of fixed temperatures at the ice 
shelf surface T s and base T B . Utilizing these boundary 
conditions, the temperature profile can be written 


T,{z) = 



+ T„ - T s exp 


w,H, 


1 - exp 



(25) 


with z — 0 at the surface and z — ~H , at the base. The 
solution expressed by (25) is shown in Fig. 2 for an ice 
shelf 1000 m thick experiencing a range of basal melting 
and freezing rates. Under melting conditions, the tem- 
perature gradient near the base of the ice shelf is in- 
creased, which increases the conductive heat flux and 
serves to counteract the melting. Conversely, freezing 
decreases the temperature gradient and the conductive 
flux, which leads to a lower freezing rate than would 
be estimated for a linear temperature profile. 

From Eq. (25) we can derive the temperature gradient 
at the ice shelf base: 


dT, 

dz 


B 


= n 


(T s - T b ) 

Hi 


(26) 


Equation (26) has the same form as Eq. (21 ) apart from 
the factor 
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n = 



(27) 


which depends on the Peclet number 


Y = 



(28) 


For a melt rate of 1 myr 1 at the base of an ice shelf 
1000 m thick we obtain a Peclet number of approxi- 
mately -30, implying that the temperature gradient am- 
plification factor II can play an important role in Eq. 
(26). 

Immediately obvious is that Eq. (27) is ill-defined for 
a melt rate of zero, a problem that may be overcome 
by rewriting the right-hand side as a power series (Arf- 
ken 1970) of the form 



* Y« 

- 2 B.h, 


n\ 


(29) 


where B tl represent the sequence of Bernoulli numbers. 
The first few Bernoulli numbers have values (Abra- 
mowitz and Stegun 1972) 


B a 


1 , 


b 2 



30’ 


* 6 = 


42’ 


(30) 


the odd labeled Bernoulli numbers, £ 2 „ + , for n = 1, 2, 
3, . . . , being identically zero. The series (29) is valid 
for Y < 27r and is therefore appropriate for small values 
of the Peclet number (i.e M a near-zero melt/freeze rate). 
In particular when Y = 0, the expression in (29) equals 
B 0 and Eq. (26) is then identical to Eq. (21), the purely 
diffusive solution with no advection. We also note that 
for sea ice Y < 1 typically, so the purely diffusive 
solution is a good approximation. 

Figure 3 illustrates the behavior of II as a function 
of the Peclet number Under conditions of moderate to 
high freezing II is close to zero, while for moderate to 
high melting it is very close to the absolute value of 
the Peclet number. This suggests that a possible sim- 
plification of Eq. (27) is 



for melting case when w B > 0 
for freezing case when w B < 0. 


(31) 


This approximation (also shown in Fig. 3) has the ad- 
vantage of linearizing Eq. (2), and hence simplifying 
the solution of the ice-ocean boundary equations. Such 
an approximation was introduced by N0st and Foldvik 
(1994) and has been adopted in the model of Hellmer 
et al. (1999). 


4) More complex models 

The next stage of complexity in the heat transport 
problem would be to introduce a vertical velocity that 


Equivalent melt rate (m yr"*) 

0.72 0.54 0.36 0.18 0 -0.18 -0.36 



Fig. 3. Temperature gradient amplification factor as a function of 
Peclet number. The dotted tine shows the approximation given in Eq. 
(31). The upper axis scale indicates equivalent melt rates for an ice 
shelf 1000 m thick. 


varies linearly from the surface to the base of the ice 
shelf. This would lead to the classical solution for the 
temperature profile in an ice column, first introduced by 
Robin (1955) for ice sheets. However, the solution in- 
volves either error functions or Dawson’s integrals, and 
we have little hope of recovering a linear version of Eq. 
(2). To use such a solution, the basal temperature gra- 
dient would have to be calculated as a separate problem 
and the result introduced directly into Eq. (6). The same 
applies to models of greater sophistication that could 
include the horizontal advection of heat. Such a model 
was used by Jenkins (1991) for a specific region of 
Ronne Ice Shelf where measurements of ice flow and 
surface temperature made the computation of a steady- 
state temperature distribution within the ice shelf pos- 
sible. 

3* Comparison of model results 

We begin by discussing the behavior of the three- 
equation formulations, which are the most widely used 
in the literature on ice shelf-ocean interactions, and we 
will initially treat the ice shelf as a perfect insulator. 
With no heat conduction into the ice shelf, the only 
external parameters that enter into Eqs. (1), (2), and (5) 
are the friction velocity and the thermal driving F*, 
defined as 

F* = F„ - (aS M + b + cp B ). (32) 

The difference between the various formulations lies in 
the specification of the exchange velocities for heat and 
salt, which are illustrated in Fig. 4. The two formulations 
that include an explicit parameterization of the boundary 
layer yield heat transfer coefficients that are approxi- 
mately linear functions of friction velocity, with values 
that generally differ by no more than about 10% (Fig. 
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V 






V 



%??£gzrx - - - - 

m yr-'. (b) Ratio of salt to heat transfer coefficients for the same formulations. 


4a). The reason for this close agreement is the domi- 
nance of the molecular term, which is the same in all 
cases, in the denominators of Eqs. (11) and (14). The 
two models having constant exchange velocities use val- 
ues of y T that are consistent with a friction velocity of 
about 0.01 m s'\ which corresponds to a mixed layer 
velocity of about 0.2 ms' 1 . Although this is higher than 
most estimates of currents associated with the thermo- 
haline circulation, it is consistent with rms currents as- 
sociated with strong tidal flow (Makinson and Nicholls 
1999). That the ratio of y , to y T is significantly less 
than one in all cases, and is relatively constant in most, 
(Fig. 4b) is evidence of the importance of molecular 
processes. Heat and salt transfer within the turbulent 
part of the boundary layer has little impact on the size 
of the coefficients except in the formulation that con- 
siders the influence of gravitational stability, and even 
in this formulation the effect is only significant under 
conditions of high thermal driving and low friction ve- 

locity. . 

A model’s response to thermal driving is determined 
by the magnitudes of both the heat and the salt transfer 
coefficients. Figures 5a— c show the influence of chang- 
ing the size of both y s and y T while keeping their ratio 
constant. Temperature and salinity differences across the 
boundary layer are set entirely by the thermal driving, 
with computed melt rates then responding linearly to 
variations in y 5 and y T , as heat and salt are transported 
with varying ease across the boundary layer. This re- 


sponse is similar to what we would anticipate for a two- 
equation formulation. However, the salinity difference 
across the boundary layer (Fig. 5c) means that the cor- 
responding temperature difference (Fig. 5b) is always 
smaller than the thermal driving. Calculated melt/freeze 
rates are therefore lower than they would be if the 
boundary salinity were assumed to be equal to the mixed 
layer salinity. There is also a slight nonlinearity in the 
response to thermal driving that arises because the sa- 
linity at the ice-ocean interface can increase without 
limit but can only decrease by ~35 psu before becoming 
completely fresh. 

The role played by salinity diffusion in determining 
the melt/freeze rates is shown more clearly in Figs. 5d-f, 
where the effect of varying the ratio of y s to y T while 
keeping the latter constant is illustrated. As the ratio 
tends to infinity, the response becomes that of a two- 
equation formulation with zero salinity difference across 
the boundary layer (Fig. 50, a temperature difference 
equal to the thermal driving (Fig. 5e), and a melt/freeze 
rate that is directly proportional to T* (Fig. 5d). De- 
creasing the salinity diffusivity increases the salinity 
difference across the boundary layer, which decreases 
the temperature difference and with it the melt/freeze 
rate while the nonlinearity in the response to 7* be- 
comes more pronounced. With negative thermal driving, 
the salinity at the boundary can grow until the freezing 
point depression balances the thermal driving, yielding 
a temperature difference of zero across the boundary 
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Fig. 5. Response of a three-equation formulation for thermal driving of +2°, 4- 1°, 0°, —1°, and -2°C plotted against 
the magnitude of the thermal exchange velocity y T , (a) melt rate, (b) temperature difference across the boundary layer, 
(c) salinity difference across the boundary layer. The ratio of the salinity exchange velocity to y T is kept at 0.04. 

Panels (d), (e), and (f) show the response of the same variables to changes in the ratio of y 5 to y,, while the latter is 
kept constant at 1.0 X 10 4 m s The dotted lines indicate ratios typical of formulations with constant coefficients 
(0.005) and those based on boundary layer parameterizations (0.04). 


layer and hence no freezing, but, given sufficient pos- 
itive thermal driving, melting can proceed even if the 
water at the boundary is completely fresh. The two val- 
ues of y s fy T found in the literature span a region of high 
sensitivity. The lower ratio gives rise to melting and 
freezing rates that are not only smaller (for the same 
heat transfer coefficient), but also show a more nonlinear 
response to thermal driving. 

The effects discussed above are illustrated quantita- 
tively in Fig. 6, which shows the melt/freeze rates com- 
puted by each of the models for a broad range of thermal 
driving. The two boundary layer parameterizations give 
rather similar results, suggesting that the precise form 
of r Turb is not critical. This is convenient, as the intro- 
duction of the stability parameter into Eq. (15) involves 
a computationally expensive iteration to derive the melt 
rates and exchange coefficients. Simply setting the sta- 
bility parameter to unity does not have a large impact 
on melt rates computed with this model, except at very 
low friction velocity (Fig. 7a). For thermal driving less 
than 0.5°C (i.e., values commonly found in nature) and 
a friction velocity greater that 0.001 (corresponding to 
a velocity of about 0.02 ms 1 ) differences between melt 


rates computed with and without the stability parameter 
differ by less than 10%. 

A possible refinement to the models discussed above 
would be the introduction of a conductive heat flux into 
the ice shelf. The influences of purely diffusive and 
constant _v^icaI Mvection/diffusion models are illus- 
trated in Figs. 7b,c. The linear temperature profile as- 
sumed by the purely diffusive model causes a net shift 
toward freezing so that a positive thermal driving is 
required for zero melting. Only in the region of the melt/ 
freeze transition, where the rates are very small, does 
this approximation have a noticeable impact on model 
results (Fig. 7b). The model with constant vertical heat 
advection in the ice shelf has no effect unless the mixed 
layer is warmer than the freezing point. It then reduces 
the computed melt rates by about 10% (Fig. 7c). 

McPhee (1992) and McPhee et al. (1999) demon- 
strated that direct measurements of the turbulent heat 
flux in the boundary layer beneath drifting sea ice could 
be well fitted with a two-equation model. This might 
be anticipated from the approximately linear response 
of the three-equation models, particularly at the mod- 
erate levels of thermal driving that are of most practical 
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Fig 6 Melt/freeze rates as a function of thermal driving calculated using exchange velocities given by 
Eos. ( 1 4)— ( 18) (solid lines), by Jenkins (1991) (dotted lines), by Hellmer and 0, bers (1989) (dashed line), 
by Scheduikat and Olbers (1990) (dot-dashed line), and by Determann and Gerdes (1994) (solid line with 
dots). In the Jenkins and Hellmer and Olbers cases, three curves are shown and labeled for friction velocities 
of 0 0 1.0, and 2.0 cm s - '. Panel (b) shows an enlargement of the boxed area in panel (a). 


interest (Fig. 6b). However, the key to a consistent two- 
equation formulation lies in the choice of an effective 
exchange velocity, which accounts approximately for 
the fact that the finite salinity diffusivity supports a 
salinity difference across the boundary layer and hence 
reduces the thermal forcing. Using Eq. (1) we can re- 
write Eq. (32) as 

T* = (T„ - T„)~ a (S M - S„) (33) 

from which we can express the ratio of thermal driving 
(7*) to thermal forcing (7 M - T B ) as 


0 = 1 - a 


(S M - S B ) 


(34) 


{T„ ~ T B ) 

From Eqs. (2) and (5), with all fluxes into the ice shelf 
ignored, Eq. (34) can be written 

L f y s 


0 = 1 - 


(35) 


Provided this factor is approximately constant, a two- 
equation formulation with an effective transfer coeffi- 
cient of y T /& should yield reasonable melt rates. Taking 
S H = 34.5 psu gives 0 = 1.6 for y s ly T = 0.04, typical 
for the models of Jenkins (1991) and McPhee et al. 
(1987), and 0 = 5.7 for y s /y r = 0.005, as in the models 
of Hellmer and Olbers (1989) and Scheduikat and 01- 
bers (1990). 

Figure 7d illustrates the differences between melting/ 
freezing rates calculated with the model that includes 


the stability parameter and those derived from an equiv- 
alent two-equation formulation. We find large differ- 
ences for high thermal driving, particularly at low values 
of the friction velocity, most of which are a result of 
ignoring the effect of stability (Fig. 7a). At higher fric- 
tion velocity, the linear response of the two-equation 
formulation means that melting rates tend to be under- 
estimated and freezing rates overestimated compared 
with the results of the full three-equation model. How- 
ever, for conditions frequently encountered in nature 
( 17*1 < 0.5°C, w* > 0.001) differences between the two- 
and three-equation formulations are typically less than 
10 %. 

In Fig. 8 we compare the effective transfer coetti- 
cients for each of the formulations discussed above, with 
that derived by McPhee et al. (1999) from observations 
in the turbulent boundary layer beneath sea ice. The 
models of McPhee et al. (1987) and Jenkins (1991) re- 
produce the observed dependency on friction velocity, 
but overestimate turbulent transfer by about 15% and 
30%, respectively. The constant transfer coefficients 
used by Hellmer and Olbers (1989) and Scheduikat and 
Olbers (1990) are consistent with currents of 0.06-0.08 
m s“ 5 , the right order of magnitude for the thermohaline- 
forced circulation beneath ice shelves. The effective 
transfer coefficient of Determann and Gerdes (1994) is 
much higher, corresponding to a velocity of about 0.35 
m s', and is therefore appropriate for a cavity subject 
to vigorous tidal mixing. The melt/freeze rates produced 
by this latter model are shown in Fig. 6 for comparison. 
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Fig. 7. Ratio of the melt/freeze rate r derived from different formulations to that calculated 
using the standard three-equation Formulation with exchange velocities set according to Eqs. (14)- 
(18) and no heat conduction into the ice shelf. In each panel, solid lines indicate results obtained 
with a friction velocity of 1 cm s ', dashed lines with 0. 1 cm s ', and dotted lines with 0.01 cm 
s _l . The different formulations are (a) exchange velocities derived from Eqs. (14)-(18) with the 
stability parameter set to l, (b) the standard model with vertical heat diffusion in the ice shelf, 
(c) the standard model with vertical advection and diffusion of heat in the ice shelf, and (d) an 
equivalent two-equation formulation. Note that the vertical scale differs between panels. 
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Fig. 8. Effective transfer coefficients for the formulations intro- 
duced by Hellmer and Olbers (1989) (dashed line labeled H + O), by 
Scheduikat and Olbers (1990) (dotted line labeled S + O), by Jenkins 
(1991) (solid line labeled J), in this paper (unlabeled, dashed line), 
and by Determann and Gerdes (1994) (solid line labeled D + G). The 
dotted line labeled MKM illustrates the effective transfer coefficient 
derived by McPhee et al. (1999) based on observations. 


4. Computed buoyancy fluxes 

Buoyancy fluxes associated with melting and freezing 
represent the primary forcing on the ocean beneath an 
ice shelf. Here we analyze how the differing formula- 
tions of the ice-ocean interaction discussed above in- 
fluence the forcing imparted to an ocean model. The 
rate of change of mixed layer buoyancy can be written 

- IS} <36 » 

where and /3 r are the salinity contraction and thermal 
expansion coefficients, respectively. The term outside 
the brackets is directly related to the freshwater flux 
associated with melting: 

S M = ~ (37) 

where H M denotes the thickness of the mixed layer. The 
term in brackets can be rewritten 

(3j M = (BT/dS)„ 

nx (dT/dsi; 

Melting of ice into the mixed layer causes its T/S prop- 
erties to evolve along a straight line, the gradient of 
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which has been derived by Gade (1979), Greisman 
(1979), and N0st and Foldvik (1994), and appears as 
the numerator on the right-hand side of Eq. (38). The 
denominator is the isopycnal slope in TfS space. Of 
relevance to the discussion here is that (dTldS) M depends 
on the assumptions made about the heat flux into the 
ice shelf (N0st and Foldvik 1994). In particular, the 
advection/diffusion model yields a ( dT/dS) M of 2.8 (for 
a surface temperature of — 25°C), while the commonly 
used model of a perfectly insulating ice shelf gives a 
value that is lower by 0.33. The overall error in buoy- 
ancy forcing is 

*»• - m + (39> 

from which we obtain 

8B^ = 8w, 8(8TfdS)„ 

K w B (dTldS) a - (dT/dS) M 

As expected, we find that any model that gives an error 
in the melt rate contributes the same percentage error 
to the buoyancy forcing. However, if the melt rate is 
misdiagnosed because of an error in the estimate of the 
heat conducted into the ice shelf, there is an additional 
error in the forcing. Figure 7c shows a 10% difference 
between the melt rates derived from a model having 
constant vertical advection in the ice shelf and those 
derived with a model that treats the ice as a perfect 
insulator. The additional error in buoyancy forcing, aris- 
ing from the second term on the right-hand side of Eq. 
(40), is small (—1%) near the surface, where the iso- 
pycnal s are steep, but rises to about 5% at a depth of 
2000 m, which is reached beneath the thickest ice 
shelves. 

5. Summary and conclusions 

The main objective of this study has been the pre- 
sentation of an hierarchy of models describing the ther- 
modynamic interaction between the base of an ice shelf 
and the underlying ocean waters. We have reviewed the 
various models that have been used in the literature on 
ice shelf-ocean interactions and have introduced a pa- 
rameterization of turbulent transfer in the oceanic 
boundary layer, based on the work of McPhee et al. 
(1987), that considers the impact of gravitational sta- 
bility. We have investigated the behavior of all the mod- 
els and analyzed their performance in the light of recent 
studies of the turbulent boundary layer beneath sea ice 
(McPhee 1992; McPhee et al. 1999). A key finding of 
the latter authors is that turbulent transfer is apparently 
independent of the roughness of the ice-ocean interface, 
a fact that gives us confidence in extrapolating their 
findings to an ice shelf base of unknown roughness. 

Most of the models in the literature conform to the 
three-equation formulation, although the choice of ther- 
mal and salinity transfer coefficients varies. We have 


shown that the behavior of these models can be ap- 
proximated by an equivalent two-equation formulation, 
at least for moderate thermal driving. The nonlinearity 
in the response of the three-equation models, a feature 
that does not appear in the simpler formulation, only 
becomes apparent at high thermal driving. In nature, 
supercooling can always be damped by the formation 
of frazil ice within the water column (Jenkins and Bom- 
bosch 1995), making model behavior at negative ther- 
mal driving greater than ~0.1°C of theoretical rather 
than practical interest. Conditions of high positive ther- 
mal driving are unlikely to be encountered beneath the 
ice shelves of the Ross and Weddell Seas, but mea- 
surements from beneath George VI Ice Shelf, in the 
Bellingshausen Sea, show water more than 1°C above 
freezing within a few meters of the ice shelf base (K. 
W. Nicholls 1998, personal communication). A two- 
equation formulation may be inappropriate under these 
conditions, which may require the use of the full model 
of McPhee et al. (1987). 

We have discussed various parameterizations of the 
heat flux into the ice shelf. The only one that manages 
to capture any of the nonlinearity of the typical ice shelf 
temperature profile is that which assumes constant ver- 
tical advection of ice. Applying this parameterization 
reduces melting by about 10% but reduces the buoyancy 
forcing on the ocean by up to 1 5%, the additional change 
being the result of the heat loss to the overlying ice. 
The overall effect is comparable to the differences in 
buoyancy forcing associated with the various choices 
of transfer coefficients used in the three-equation mod- 
els. In nature, the temperature distribution within an ice 
shelf is determined by the history of melting and freez- 
ing that each ice column has experienced, but to intro- 
duce this thermal memory of past conditions would re- 
quire a rather sophisticated dynamic/thermodynamic 
model of the ice shelf. 

Our most important results are summarized in Figs. 
6 and 8, which illustrate the behavior of the models 
used to date in the literature on ice shelf-ocean inter- 
actions. The three-equation models of Hellmer and Ol- 
bers (1989) and Scheduikat and Olbers (1990) have ef- 
fective transfer coefficients that are only one-fifth the 
size of that used by Determann and Gerdes (1994). In 
the models of Jenkins (1991) and Jenkins and Bombosch 
(1995), typical friction velocities lie in the range 0.002- 
0.008 ms”’, yielding effective transfer coefficients 
ranging in magnitude from that of Hellmer and Olbers 
(1989) up to half that used by Determann and Gerdes 
(1994). Whether it is better to use constant coefficients, 
ones based on assumed tidal velocities, or ones based 
on computed thermohaline velocities is an open ques- 
tion. However, when comparing model output, it is im- 
portant to realize that the differing parameterizations of 
the ice— ocean interaction yield melting rates and hence 
buoyancy fluxes, that vary over a factor of 5 for the 
same thermal driving. Our recommendation is that for- 
mulations used in future, whether two-equation or three- 
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equation, should aim to reproduce the behavior ob- 
served by McPhee (1992) and McPhee et al. (1999), at 
least until such time as measurements of turbulent heat 
flux beneath ice shelves are available. 
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Introduction 

The seawater-filled cavity beneath Filchner-Ronne Ice Shelf (FRIS) is a region of major 
importance for both the Antarctic Ice Sheet and the Southern Ocean. High Salinity Shelf Water 
(HSSW) flowing into the cavity melts ice from the underside of the ice shelf and in the process 
is transformed into Ice Shelf Water (ISW), a precursor of the coldest and densest form of 
Antarctic Bottom Water. However, direct observations of conditions in the cavity are restricted 
to a few isolated sites and there remains considerable uncertainty over the strength of the 
exchange between the cavity and the open ocean to the north. The production of HSSW over 
the continental shelf of the southern Weddell Sea varies both spatially and temporally, but the 
consequences of this variability for the ice shelf-ocean system are presently unknown. 

A study of water properties measured near the ice front led N0st and Foldvik ( 1 994) to the 
conclusion that ISW in the Filchner Depression is formed from the HSSW found in the Ronne 
Depression. This indicated a west-to-east circulation beneath the ice shelf, but Nicholls et al. 
(1997) suggested that most of the inflowing water reaching the deeper parts of the Ronne 
Depression is unable to cross the relatively shallow seabed to the south of Korff and Henry ice 
rises and must remain trapped within the depression. More recent observations reported by 
Nicholls et al. (this volume) suggest strong west-to-east flow around the southern tip of Berkner 
Island, but the western source region is not clear. Currentmeter measurements discussed by 
0sterhus at FRISP ‘99 suggest that some may come from the Berkner Bank. Results from a 
two-dimensional model of thermohaline circulation beneath FRIS indicated the possibility of 
flow around Berkner Island (Hellmer and Olbers, 1991), but more recent three-dimensional 
studies suggest only minor west-to-east exchange, most of the flow being concentrated within 
gyres that are strongly constrained by the sub-ice topography (Grosfeld and Gerdes, 1998; 
Gerdes et al., 1999). 

Here we describe the first results obtained with an isopycnic coordinate ocean model. Our aim 
is to investigate the response of the cavity beneath Filchner-Ronne Ice Shelf to the formation of 
HSSW on the continental shelf to the north. The model domain therefore includes the whole of 
the continental shelf and is forced only by the restoration of surface temperature and salinity to 
a prescribed seasonal cycle over the open ocean portion of the domain. 
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The model 

We have adapted the Miami Isopycnic Coordinate Ocean Model (MICOM; Bleck et al., 1992) 
so that it can be used for sub-ice-shelf investigations. The key feature of this model is that the 
vertical discretisation of the domain is based on density. The ocean is divided into a number of 
layers of constant-density water, and the model computes the thickness of these layers at each 
point of the horizontal grid. The advantage of this approach is that any artificial, numerical 
diffusion that is present is purely isopycnal. Diapycnal diffusion can be carefully controlled and 
even eliminated altogether, if desired. For the model to accept buoyancy forcing, the upper layer 
must have a freely-evolving density, so MICOM uses an embedded mixed layer at the surface. 
The entrainment/detrainment algorithm allows water to be cycled from one isopycnic layer, 
through the mixed layer, where its properties are modified by the surface boundary conditions, 
into another isopycnic layer. Unlike the isopycnic layers that can appear and disappear as 
required, the surface mixed layer is always present at all grid points, including those beneath the 
ice shelf. Most of the modifications we have made to the code have been connected with the 
mixed layer routines. 

The model domain (Figure 1) covers 12° of latitude and 60° of longitude and is divided into 
square grid cells having dimensions equivalent to 0.75° of longitude. This gives a horizontal 
resolution ranging from 1 1.5 km along the southern boundary, to 28.5 km along the northern 
boundary. Ice thickness and seabed depth south of 75°S are taken from Vaughan et al. (1995), 
while to the north seabed depths come from the ET0P05 database. Initial conditions are taken 
from the AWI Hydrographic Atlas (Olbers et al., 1992), with the cavity being filled by 
extrapolation. The northern and eastern boundaries are solid walls, where scalar variables are 
restored to Altas values. The only external forcing on the domain is the restoration of 
temperature and salinity in the mixed layer to the north of the ice front. Winter and summer 
salinity fields are shown in Figure 2. At each point the variation between the two is a simple 



Longitude 


Figure 1 : Water column thickness (m) over the southern Weddell 
Sea domain. 
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Figure 2a: Mid-summer restoring field for Figure 2b: Mid-winter restoring field for sea 
sea surface salinity (psu). surface salinity (psu). 


sinusoid. Temperature is restored to the surface freezing point during most of the year, but is 
ramped up to -1°C in December and January. 


Preliminary results 

The model has been run out over ten years with the forcing described above. While it would 
appear that some form of steady state has been reached at this stage (Figure 3), it should be noted 
that the south and west of the cavity are rather poorly ventilated, and may still by influenced by 
the initial conditions. In general, the computed melt/freeze rates accord reasonably well with 



Figure 3: Melt/freeze rates (cm yr' 1 ) averaged 
over the ice shelf base. The upper line shows the 
mean of all grid points experiencing melting and 
the lower line the mean of all points experiencing 
freezing. The centre line is the overall average of 
all points in the sub-ice-shelf domain. 
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those deduced from observation. Melting is focussed near the grounding lines of the major ice 
streams, as we would anticipate, but outside of these regions mass exchange between ice and 
ocean appears to be too weak. This is likely an artefact of the mixed layer 
entrainment/detrainment routine, which represents the only mechanism by which heat can be 
brought to the ice shelf base. The peak melt rates in mid-summer (Figure 3) reflect the influence 
of the warm surface layer on the regions near the ice front. The secondary peak in late winter, 
at a time when there is an associated peak in freezing, corresponds to the time of maximum 
external forcing on the cavity circulation. The overall average is a net melt rate of 20 cm yr 

The impact of the external forcing, which is strongest in late winter when convection over the 
open continental shelf reaches its deepest, can clearly be seen in Figure 4. This figure also 
illustrates that the depth-averaged flow produced by MICOM differs markedly from that reported 
by Grosfeld and Gerdes (1998) and Gerdes et al. (1999). We see generally weak and variable 
depth-averaged velocities, which do not form into the strong, coherent gyres described by the 
above authors. This is because the direct impact of the buoyancy forcing is confined to the mixed 
layer. The lower layers move only in response to changes in the mixed layer, and the induced 
motion is often in a direction that differs from that of the mixed layer itself. Bearing this in mind, 
we should be cautious in interpreting Figure 4, which only illustrates one aspect of the circulation. 
However, it is clear that the cavity is ventilated primarily via the Ronne and Filchner depressions, 
and that waters reach the remotest parts of the cavity only indirectly and intermittently. 

The biggest seasonal changes in flow occur over the Filchner Depression, and these are a direct 
result of the surface forcing. In summer, a cyclonic gyre within the depression penetrates only 
weakly into the cavity, but in winter, convection on the Berkner Bank produces a mixed layer that 
extends to 800m on the western flank of the depression and induces a strong west-to-east density 
gradient (Figure 5b). The cyclonic gyre is displaced to the east and strong southward flow 
develops in the west of the depression. After mid-winter the surface buoyancy forcing weakens 
and the vertical isopycnals subside. Dense waters drain from the Berkner Bank into the 
depression, turning to the north, while mid-depth waters move to the west and turn beneath the 
ice shelf. By mid-summer most of the potential energy has been removed from the water column, 
and we are left with the domed isopycnals, associated with the cyclonic gyre that has now re- 
established itself (Figure 5a). 



Figure 5a: Late summer potential density (kg Figure 5b: Late winter potential density (kg 

m' 3 - 1000) distribution at 77°S in the Filchner m' 3 - 1000) distribution at 77°S in the Filchner 

Depression. Depression. 
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The picture of circulation described above is sensitive to the prescribed surface restoring fields. 
In earlier simulations we used a surface salinity field that differed significantly from the one 
shown in Figure 2. Lower salinities in the east meant that convective overturning never occurred 
on the Berkner Bank, and the dominant inflow to the cavity originated in the Ronne Depression. 
The inflowing waters turned to die east, followed the c oast of Berkner Island and finally emerged 
in the Filchner Depression. The area to the south of Korff and Henry ice rises was also poorly 
ventilated. The circulation described earlier should therefore not be interpreted as a definitive 
statement on the current circulation beneath Filchner-Ronne Ice Shelf, but more as the first in a 
number of sensitivity studies into the response of the cavity to the spatial distribution of HSSW 
production on the open continental shelf. 


Conclusions 

Preliminary runs of MICOM on a domain resembling the real topography of the southern Weddell 
Sea illustrate the sensitivity of the circulation beneath Filchner-Ronne Ice Shelf to the seasonal 
formation of HSSW to the north of the ice front. Although the circulation throughout the cavity 
responds to the seasonal forcing, the remotest parts of the cavity remain poorly ventilated. The 
possibility exists for waters to reach the interior from either the Ronne or the Filchner Depression, 
depending on the location of convective activity on the open continental shelf. This conclusion 
is similar to that reached by Hellmer and Olbers (1991), but the reasoning is subtly different. The 
model of Hellmer and Olbers (1991) responded directly to the density of the waters within the 
two depressions. In the results described in this paper, the densest HSSW is formed in the Ronne 
Depression, but stronger inflow develops in the Filchner Depression because of the more 
pronounced zonal density gradients generated there and the greater depth range over which these 
gradients persist. Hence, there is more potential energy available to drive shear flow 
perpendicular to the ice front. These results suggest that changes in the production of HSSW, 
such as those described by Npst and 0sterhus (1998), could have a major impact on the 
circulation within the cavity. We plan to investigate this by studying the response of our model 
to reduced restoring salinities over the Berkner Bank. 
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ABSTRACT 

Much of the Antarctic coastline comprises large, floating ice shelves, beneath which waters from the open 
ocean circulate. The interaction of the seawater with the base of these ice shelves has a bearing both on the 
rate at which Antarctic Bottom Water is formed and on the mass balance of the ice sheet. An isopycnic coordinate 
ocean general circulation model has been modified so as to allow the incorporation of a floating ice shelf as an 
upper boundary to the model domain. The modified code admits the introduction of an arbitrary surface pressure 
field and includes new algorithms for the diagnosis of entrainment into, and detrainment from, the surface mixed 
layer. Special care is needed in handling the cases where the mixed layer, and isopycnic interior layers, interact 
with surface and basal topography. The modified model is described in detail and then applied to an idealized 
ice shelf-ocean geometry. Simple tests with zero surface buoyancy forcing indicate that the introduction of the 
static surface pressure induces an insignificant motion in the underlying water. With nonzero surface buoyancy 
forcing the model produces a cyclonic circulation beneath the ice shelf. Outflow along the ice shelf base, driven 
by melting of the thickest ice, is balanced by deep inflow. The abrupt change in water column thickness at the 
ice shelf front does not form a barrier to buoyancy -driven circulation across the front. 


1. Introduction 

Eustacy is affected by a variety of factors including 
the storage of land surface and ground water, the thermal 
expansion of the oceans, and the mass balance of gla- 
ciers and ice sheets. A report by the Intergovernmental 
Panel on Climate Change (Houghton et al. 1996) has 
singled out the ice sheets as representing the greatest 
uncertainty in accounting for past, and predicting future, 
changes in sea level. Ice sheet mass balance is a com- 
bination of surface accumulation and ablation, calving 
of icebergs from the fronts of ice shelves and tidewater 
glaciers, and melting and freezing that occurs at the ice- 
ocean interface. The latter component of the mass bud- 
get has no direct impact on sea level, as the vast majority 
of the melting ice is already freely floating in the ocean, 
but the ice shelves may play a role in regulating the 
discharge from some grounded portions of the ice sheet. 

Aside from this indirect impact on global sea level, 
basal melting also has an important influence on the 
properties of the waters that circulate beneath an ice 
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shelf. Outside the subice cavity the only heat sink avail- 
able to the ocean is the atmosphere, so seawater cannot 
be cooled below the surface freezing point of about 
~ 1.9°C. The base of an ice shelf represents a heat sink 
at depths approaching 2000 m where the freezing point 
may be as low as -3.4°C (Millero 1978), so the waters 
in contact with the ice shelf may be up to T5°C below 
the freezing point at atmospheric pressure. Such poten- 
tially supercooled water is referred to as ice shelf water 
(ISW). In the Weddell Sea, ISW contributes to the for- 
mation of Antarctic Bottom Water (Foldvik and Kvinge 
1974), and the exceptionally cold temperatures attained 
by ISW enhance the density of the product water, 
through the thermobaric effect. A further consequence 
of the generation of potentially supercooled water is that 
if it rises toward the surface it may become supercooled 
with respect to the in situ freezing point. The result is 
ice production deep in the water column (Foldvik and 
Kvinge 1974). Where this happens beneath the ice shelf, 
the ice accumulates on the base as marine ice (Bom- 
bosch and Jenkins 1995), while the freezing of ISW that 
exits the subice cavity may contribute to the sea-ice 
budget (Bombosch 1998). 

The interaction between ice sheets and oceans is thus 
an important element of the climate system, and in re- 
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cent years numerical models have been used to evaluate 
the key processes operating in the subice cavity (Wil- 
liams et a!. 1998). The approach taken has been to use 
ocean models of varying sophistication and apply upper 
boundary conditions derived from a thermodynamic 
model of the ice shelf-ocean interaction (Holland and 
Jenkins 1999). The calculated heat and freshwater fluxes 
represent the only surface forcing on the sub-ice shelf 
water column, since the overlying ice shelf isolates the 
ocean from wind, evaporation and precipitation, and the 
direct effects of atmospheric temperature variations. Dy- 
namic models of the ice shelf itself have not been in- 
cluded to date, and the implicit assumption has been 
made that any ice shelf thickness changes resulting from 
melting or freezing are offset exactly by the flow of the 
ice shelf itself. The disparity of timescales between the 
slowly evolving ice shelf and the rapidly ventilated wa- 
ters beneath provides some justification for this mod- 
eling approach. 

In this paper we focus on the mathematical descrip- 
tion of an isopycnic coordinate sub-ice shelf ocean 
model, and in particular the treatment and parameteri- 
zation of the mixed layer processes. Observations show 
the presence of a mixed layer in the water column be- 
neath an ice shelf (Nicholls and Makinson 1998) with 
water properties relatively well mixed near the ice shelf 
base and sometimes over a significant fraction of the 
water column. Section 2 presents the conservation equa- 
tions for the isopycnic layers of the model and describes 
how the equations differ for the mixed layer. The nec- 
essary modifications to the model required for it to func- 
tion in the oceanographic regime beneath an ice shelf 
are then presented. The main issues are that of (i) the 
incorporation of an arbitrary surface pressure field, (ii) 
the processes of mixed layer entrainment and detrain- 
ment, and (iii) the interaction of layers with topography. 
We also briefly present some additional changes we have 
made to the code, associated with internal friction, that 
are not directly related to the ice shelf problem. Finally 
the modified model is applied to an idealized ice shelf- 
ocean geometry and the results are discussed. 

2. Model description 

As a longer-term science objective a coupled numer- 
ical model is being developed suitable for investigation 
of oceanographic processes in polar regions, referred to 
as the Polar Ocean Land Atmosphere Ice Regional (PO- 
LAIR) modeling system. At present this system consists 
of a viscous-sublayer parameterization describing the 
thermodynamic interaction between the ice shelf base 
and the ocean surface (Holland and Jenkins 1999) cou- 
pled to an isopycnic ocean general circulation model 
(Bleck et al. 1992). That ocean component, the Miami 
Isopycnic Coordinate Ocean Model (MICOM), already 
contains an embedded mixed layer model based on a 
turbulent kinetic energy (TKE) budget. In the next stag- 
es of development, a sea-ice, atmospheric boundary lay- 


er, and dynamical ice shelf model will be added. We 
now briefly describe the ocean model and its mixed layer 
submodel prior to elaborating upon the modifications 
that are addressed in the sections following. 


a. Isopycnic interior model 

The primitive equation isopycnic coordinate ocean 
model used (MICOM) has four prognostic equations: 
one equation for the horizontal velocity vector, a layer 
thickness tendency equation, and two conservative 
equations for the buoyancy related variables. Since not 
all layers are necessarily constant density (e.g., the 
mixed layer), we first present the equations for a gen- 
eralized vertical coordinate £ (Bleck and Boudra 1981; 
Bleck 1998): 


dv 

dt f 
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In these relations, v = ( u , u) is the horizontal velocity 
vector, p the pressure, C represents either one of the 
model’s dynamically arrive tracers (i.e., potential tem- 
perature 0 or salinity S) as well as any passive tracers, 
£ = dvfdx^ — du/dy 4 the relative vorticity, a 9 the specific 
volume, <f> = gz the geopotential, g the gravitational 
acceleration, / the Coriolis parameter (spatially vary- 
ing), k the vertical unit vector, k v the isopycnic eddy 
viscosity, K h the layer thickness diffusivity, k c the tracer 
diffusivity, t the wind- and/or bottom-drag-induced 
shear stress vector, and J c represents the sum of all 
diabatic forcings on the tracer variable C. The quantity 
(dp/d£ B#dt) represents the vertical mass flux through a 
surface £ The independent variables x, y, and t have 
their usual meanings and £ is now the independent var- 
iable in the vertical direction. 

The above equations are relationships formulated for 
variables on £ surfaces, but the ocean model is for- 
mulated for layers (i.e., the material bounded between 
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two £ surfaces). A point of nomenclature: in the discrete 
context the individual layers will be denoted by integer 
subscript k f with the topmost layer assigned index 1 and 
the bottommost layer assigned index N . Furthermore, 
in general we will use the notation of a superscript sym- 
bol (— ) to refer to a quantity at the surface interface of 
a layer and a superscript symbol (+) at the basal in- 
terface of a layer. To obtain averaged equations for the 
fluid between some upper and lower bounding 
surfaces requires multiplication of each of (2.1) through 
(2.3) by dp/d£ vertical integration of the resulting equa- 
tions over the coordinate interval A£ = and 

division by the pressure jump between the surfaces, A p 

- “ P(€~)- As a further simplification for use in 
later sections, we introduce the notation that the layer 
thickness is exactly defined in terms of the pressure 
jump as h = A p. As the velocity within a layer is ver- 
tically constant, it is convenient to also introduce the 
notation of Au~, Au + as the velocity jump across the 
upper — and lower + interfaces of a layer, respectively. 

After these transformations, the horizontal momen- 
tum, layer thickness, and; buoyancy conservation equa- 
tions for a layer become 

y 2 
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(2.5) 
C 

( 2 . 6 ) 


The topmost layer is not isopycnic, is referred to as the 
mixed layer, and has a nonzero (dpfd^ d£/dt)f repre- 
senting a large flux of mass and properties into and out 
of an otherwise almost-adiabatic ocean interior. The in- 
terior ocean is not exactly adiabatic as there exists subtle 
diapycnic fluxes (dpldg d£Jdt)l between neighboring 
layers due to the background presence of small-scale 
turbulence. A detailed treatment of these terms is pre- 
sented elsewhere (McDougall and Dewar 1998) and no 
modification of that treatment is undertaken in this 
study. 

For all interior model layers (i.e., excluding the mixed 
layer) we choose f to be surfaces of constant potential 
density p 0 , and for these layers the (dp/di; d£/dr); terms 


vanish if there are no diapycnic exchanges. For future 
reference we will define this flux of mass between layers 
as a general diapycnic flux by introducing the notation 
wf — (dpld£ d£fdt)t . On isopycnic layers the pressure 
and geopotential gradient terms in (2.4) can be replaced 
by the gradient of the Montgomery potential (Mont- 
gomery 1937): 


M = ot d p 4- gz. (2.7) 


Because the Montgomery potential is vertically invari- 
ant throughout the depth of a layer, it does not matter 
at what level within the layer it is evaluated. For adi- 
abatic isopycnic interior layers Eqs. (2.4)-(2.6) reduce 
to (Bleck 1998) 
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As these equations are derived, whether for noniso- 
pycnic or isopycnic layers, under the hydrostatic (and 
Boussinesq) approximation, the vertical component of 
the momentum balance is a diagnostic equation, which 
takes on the form (Bleck 1998) 


dM 

da & 


= P- 


( 2 . 11 ) 


b. Mixed layer model 

A characteristic problem of pure isopycnic coordinate 
formulations is an inability to accommodate buoyancy 
forcing at the ocean surface. The solution adopted in 
MICOM is to have a variable density mixed layer sitting 
on top of the pure isopycnic layers that acts as an in- 
terface between the adiabatic ocean interior and the 
buoyancy forcing from the atmosphere. In the present 
application the surface buoyancy fluxes also arise from 
thermodynamic interactions at the base of an ice shelf, 
and so the mixed layer must also exist beneath the ice 
shelf. While all interior isopycnic layers adhere to the 
basic conservation equations (2.4)-(2.6), the mixed lay- 
er is unique in that its layer thickness is also dependent 
upon a TKE balance, determined from the input of fric- 
tional stress and buoyancy at the surface, and a param- 
eterization of energy dissipation processes (Gaspar 
1988). We define the buoyancy in a layer, b, as 

b = 5^—^. (2.12) 

Po 

where the potential density is p 6 = \la B and p„ = 1 fa a 
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is a reference density based on a reference specific vol- 
ume a 0 . Then, the rate of entrainment of fluid into the 
mixed layer, wf, can be expressed by 

= 9.W& + ^ _ (2 . I3 ) 

2 Peg &P\ 2 

where A 6 is the change in mixed layer buoyancy re- 
sulting from entrainment, m g is an empirical constant, 
w* is the friction velocity, A p, is the thickness of the 
mixed layer in pressure units, B M is the buoyancy flux 
due to surface forcing, and e d is the parameterized dis- 
sipation of TKE (see Caspar 1988 for further details). 
The relation between (2.4)-(2.6) and (2.13) is mani- 
fested through the common definition of the diapycnic 
velocity, wf, which permits mass exchange between the 
mixed layer and the isopycnic interior layers by moving 
the interface representing the mixed layer base. 

A source of TKE (e.g., strong shear in the flow be- 
tween the ocean surface and the base of the ice shelf 
or alternatively brine rejection associated with marine 
ice growth) leads in general to a deepening of the mixed 
layer; the reverse being true for a sink of TKE (e.g., 
input of meltwater during a basal melting event). During 
times when the mixed layer is deepening, water is en- 
trained from the adiabatic ocean interior into the mixed 
layer. The exact amount of water entrained is diagnosed 
by balancing the supply of TKE against the increase in 
potential energy of the water column resulting from 
mixing. At times when the mixed layer is shallowing, 
water is detrained from the mixed layer into the interior 
isopycnic layer closest in density to that of the detrained 
water. In this fashion, the waters in the oceanic interior 
are ventilated by the mixed layer. 

c. Equation of state 

The equation of state we utilize to compute potential 
density p 0 , and its inverse the potential specific volume 
a d , is an approximation of the United Nations Educa- 
tional, Scientific and Cultural Organization (UNESCO; 
Millero et al. 1980) equation of state in the form of a 
polynomial that is cubic in potential temperature, qua- 
dratic in pressure, and linear in salinity (Brydon et al. 
1999). This polynomial represents a reasonable com- 
promise between accuracy and computational economy 
as compared to the much higher-order UNESCO poly- 
nomial. It also allows for direct analytical “inversion” 
of the equation of state so that, for instance, potential 
temperature can be obtained from given values of po- 
tential density, salinity, and pressure. The polynomial 
equation for the density is, expressed in sigma-theta 
terms, 

<r(0, S, Po ) = cf(/?J + c;(p 0 )6 + cf(p 0 )S 

+ cZ(p 0 )6 2 + c%{p o )S0 + cZ(p o )0 3 
+ c?(p o )S0\ (2.14) 


where the empirical coefficients C', r • • C , are known 
functions of ocean reference pressure p n (Brydon et al. 
1999). This equation of state is suitable over the range 
of temperatures, salinities, and pressures encountered in 
a sub-ice shelf cavity. 

3. Surface pressure 

a. Isostatic assumption 

A common assumption made in the construction of 
numerical models of ocean circulation is that the surface 
pressure either is spatially and temporally constant, or 
varies at most by a fraction of the overlying atmospheric 
pressure. The ocean surface then deviates only slightly 
from a surface of constant geopotential and the devia- 
tions are always small in comparison with the depth of 
the water column. However, a floating ice shelf may 
exert a pressure approaching 200 atm on the underlying 
ocean, and the associated surface topography is then of 
the same order of magnitude as the ocean bottom to- 
pography. Computing the oceanic response to changes 
in pressure of this magnitude would be a challenging 
problem, but the ice shelf thickness evolves on time- 
scales that are much longer than those of interest for 
the dynamic response of the ocean. We can therefore 
assume that the surface pressure is static and that the 
surface elevation is always close to that of the isostat- 
ically adjusted sea level. The isostatically adjusted sea 
surface then plays a role identical to that of a constant 
geopotential sea surface in conventional ocean models. 

In the case of MICOM, the isopycnic interior layers 
are bounded by interfaces that experience a spatially 
and temporally varying pressure field and may have an 
entirely arbitrary interfacial topography. The bottom in- 
terface of the bottommost layer is unique in that it has 
a fixed bottom topography and a bottom pressure that 
undergoes small, temporal variations according to the 
total weight of fluid in the overlying water column. The 
surface interface of the mixed layer is also unique in 
that it has a small, temporally and spatially variable 
surface topography and a zero surface pressure. We un- 
dertook to transform the model code to allow for an 
arbitrarily large, static surface pressure field, as might 
be imposed by a floating ice shelf, which provides the 
surface interface of the mixed layer with a large, spa- 
tially variable surface topography. The surface of the 
mixed layer continues to undergo small height pertur- 
bations about its mean topographic profile — that profile 
now being the bottom interface of an ice shelf. After 
the modifications, there still remains no temporal var- 
iation in surface pressure, whatsoever. 

We define the isostatically adjusted sea level zf by 

(3.1) 

where a x is the specific volume of the mixed layer. In 
the open-ocean part of the domain p x is considered as 
the sea surface pressure while in the sub-ice shelf cavity 
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Fig. L Side view of horizontally varying surface pressure over the 
ocean surface due to the presence of an ice shelf. The ice shelf has 
horizontally varying thickness H , and it ends abruptly at the ice shelf 
front, creating the probable situation of the mixed layer becoming 
horizontally discontinuous (i.e., a cutline) at the ice shelf front. The 
top surface of the ice shelf and of the open ocean are denoted by 
elevation and pressure z s and p s . The surface of the ocean (i.e., the 
surface of the mixed layer) has elevation and pressure denoted by z t 
and p lt while the bottom of the mixed layer is denoted by z 2 and p 2 . 
The heavy dashed line shows for reference the equilibrium elevation 
for an ocean surface in the absence of an ice shelf that would be, of 
course, a horizontally uniform geopotential surface. 


part of the domain it is the ice shelf bottom pressure. 
We estimate the surface pressure according to 


where a, and H, are the mean specific volume and thick- 
ness of the overlying ice column, respectively. A sche- 
matic of the definition of these and other relevant quan- 
tities is given in Fig. 1. The expression (3.1) is based 
on the assumption that the ice shelf cannot support ver- 
tical shear stresses between neighboring columns of ice, 
and is a good approximation over horizontal scales that 
are large compared with the ice thickness (~1 km) as 
is the case in the present application. The sea level z? 
so defined is the equilibrium sea surface position and 
is not a geopotential surface as is typically the case in 
ocean modeling. As dyna mic al motions occur in the 
ocean, the model sea level, Z\, undergoes small height 
perturbations relative to its reference value zf« 

An important consequence of the definition of the sea 
surface (3.1), whether it be that residing beneath an ice 
shelf or that occurring over the open ocean, is that the 
gradient of the Montgomery potential will correctly van- 
ish for the situation where the horizontal pressure force 
is zero and no flow is expected. Such a situation being 
when there is no horizontal gradient of the mixed layer 
specific volume a, and the sea surface elevation is equal 
to its reference value of zf . In that case VM , vanishes 
exactly as seen from combining (3.1) with the definition 
of the Montgomery potential. 


b. Operator splitting 

In common with many ocean models, MICOM em- 
ploys an operator splitting technique to solve the mo- 
mentum equations in a computationally efficient man- 
ner. The solution for the depth-averaged flow is ad- 
vanced separately from that for the more complex, 
depth-dependent component. The former is forced by 
pressure deviations at the seabed, which reflect net con- 
vergence or divergence of mass in the overlying water 
column, while the latter is forced by changes in the 
vertical density structure, which leave the overall mass 
of the water column, and hence the pressure at the sea- 
bed, unaffected. The details of the operator splitting are 
described by Bleck and Smith (1990) and Higdon and 
Bennett (1996), so only the essential modifications will 
be presented here. The total pressure p is multiplica- 
tively split into “baroclinic” and “barotropic” com- 
ponents according to 

p = p'{ 1 4- 77 ), (3-3) 

where p f represents the depth-dependent baroclinic 
component and the dimensionless variable 77 represents 
the depth-independent barotropic component as the frac- 
tional change in total mass of the water column. At the 
seabed, the baroclinic component p f b is time invariant, 
while spatial gradients in the evolving barotropic com- 
ponent p’ h rj provide the main forcing on the depth-av- 
eraged flow. In the presence of a static surface pressure, 
Eq. (3.3) becomes 

p - p x = (p’ - p.)(1 + vY ( 3 - 4 ) 

Note that the surface pressure does not decompose, be- 
cause it is unaffected by any changes (barotropic or 
baroclinic) in the ocean. The velocity is split according to 

v = v 4 - v\ (3.5) 

where v is the depth-mean flow and the deviations v' 
have the property that they integrate to zero over the 
total water column. 


c. Barotropic component 

The prognostic equation for the depth-averaged ve- 
locity is written 

— + f k x v = — VAf 4 — - , (3.6) 

dt dt 

where we have introduced and defined notation for the 
barotropic component of the Montgomery potential as 
M = a 0 {Pb “ P\)n The final term on the right-hand 
side of (3.6) represents all addition forcing terms in (2.8) 
that have not been explicitly taken into account in (3.6). 
Applying the pressure and velocity splitting to the layer 
thickness tendency equation (2.9), we obtain 
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V ~ ^ + V ■ (vA p’) + V • (v'A p') 

= V ■ KVA//), (3.7) 

where we have used the facts that A p = Ap'(l + 77 ) 
and that 77 C 1 so that 1 + 17 may be approximated by 
1 where it appears as a factor Summing this expression 
over all layers and recalling that p[ and are constants, 
not subject to time evolution or diffusion, and that the 
velocity deviations must integrate to zero, yields a pre- 
dictive equation for the barotropic pressure component: 


The reference to the geopotential function is eliminated 
by substitution of the definition of the Montgomery po- 
tential so as to express VA/f, only in terms of quantities 
explicitly calculated by the model. The final expression 
for the baroclinic component of the mixed layer Mont- 
gomery potential is then 


VA/f = a,V 




+ V 



+ M’ 2 ) - (a l p l + a 2 p 2 ) 
2 


(3.13) 


dV + V • [V(P£ ~ Pi)1 
<>1 pi ~ Pi 


where M[ and M \ are known from (3.1 1) and the in- 
(3.8) t er f ace pressure p 2 is known from an evaluation of (3.9). 


d. Baroclinic component 

The corresponding equations for the depth-dependent 
layer thickness tendencies is obtained by substituting 
(3.8) into (3.7), 


<9A p' 
dt 


V • (v'A//) 


A// 

Pi “ P 1 


v ■ WpL 


p 1 )] 


- V - (vA p), (3.9) 

and for the depth-dependent layer momentum by sub- 
tracting (3.6) from (2.8), 

J + V^ + (£ + /# x v' + # x ? 

dt Z 

= -VA/' - — (r + - t") 

A p' 

1 dv* 

+ — V ■ KAp'Vv) - — , (3.10) 

Ap at 

where AT is the baroclinic component_of the Montgom- 
ery potential, defined by M r = M — M . The baroclinic 
Montgomery potential is computed for each layer, in- 
cluding the mixed layer, by applying the isopycnic form 
of the hydrostatic relation, Eq. (2.1 1). This is achieved by 
starting at the bottom layer, which has a known reference 
Montgomery potential, and integrating vertically, as 

Ml = AC, + («* - - A/> (3.1 1) 

where again k denotes the layer index, with k = 1 cor- 
responding to the mixed layer and k > 1 corresponding 
to the isopycnic interior layers. 

Equation (3.10) cannot be directly applied to the 
mixed layer because the density varies horizontally. We 
must therefore retain the pressure and geopotential gra- 
dients as separate terms, as in (2.4). The forcing term 
equivalent to VA/' in (3.10) for the mixed layer, denoted 
VA/f, is 


VA/f — a,V 


Pi + Pi 


+ V 


<t> 1 + 4>i 


(3.12) 


4. Mixed layer processes 

a. Entrainment 

An elegant algorithm for entrainment in an isopycnic 
layer context, consistent with overall energy conser- 
vation in the water column, is presently employed in 
the MICOM code (Bleck et al. 1989). The TKE pro- 
duced via the surface stress and buoyancy extraction 
processes is used to mix denser water from the isopycnic 
interior into the mixed layer, thereby increasing the 
mixed layer thickness. The resulting increase in poten- 
tial energy of the water column exactly matches and 
thus exhausts the supply of TKE over the course of one 
model time step. The entrainment algorithm we use is 
based on identical principles, but looks algebraically 
different, as we have had to generalize the equations to 
allow for an arbitrary surface pressure field. The deri- 
vation follows closely that of Bleck et al. (1989) and a 
schematic of the relationship between the various den- 
sities and combinations of surface buoyancy fluxes for 
a typical entrainment scenario is presented in Fig. 2. 

A point of nomenclature that needs be mentioned, as 
it applies both to the entrainment procedure of the pre- 
sent section and detrainment procedure of the next sec- 
tion, is that in the overall context of the model’s mixed 
layer algorithm there always exist three mixed layer 
thicknesses: the initial mixed layer thickness /i, at the 
beginning of a time step, an intermediate mixed layer 
thickness h g derived from Eq. (2. 13) based on the phys- 
ics of Gaspar (1988) and referred to as the “Gaspar” 
thickness, and the final mixed layer thickness h m at the 
end of a time step. 

We need to evaluate the potential energy (PE) of the 
water column between the two pressure interfaces rep- 
resenting the sea surface /?, and the unknown depth to 
which mixing and entrainment occurs p m . This will be 
done for both the mixed and unmixed states, with the 
aim of finding the value of p m for which the change in 
PE equals the supply of TKE. Introducing the gravi- 
tational acceleration as an explicit factor, the PE per unit 
horizontal area can be expressed in terms of the Mont- 
gomery potential: 
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Fig. 2. Schematic of the entrainment of isopycnic-layer waters into the mixed layer due to 
either free convection as promoted by an unstable vertical density profile or forced convection 
due to either surface frictional stress u? or an upward-directed surface buoyancy flux B m (i.e., 
buoyancy loss). The leftmost panel shows the initial state of the system with the mixed layer 
density, temperature, salinity, and thickness denoted <r,, 0 lf and h lt respectively. The com- 
plementary properties of the isopycnic layer below the mixed layer are indicated by the k 
subscript. In the presence of an upward surface buoyancy flux B m , the mixed layer properties 
are transformed into so-called Gaspar layer properties as shown in the middle panel, and denoted 
by subscript g. This is an intermediate state of the system in that the buoyancy fluxes are 
applied only over the initial mixed layer thickness, i.e., h x = /j, . This intermediate state includes 
the possibility that the mixed layer waters have become convectively unstable with respect to 
the isopycnic layer waters below. The rightmost panel shows the final state after entrainment 
has occurred with the new, deepened mixed layer overlying the partially depleted isopycnic 
layer. The final mixed layer properties are denoted with m subscripts and the final isopycnic 
layer properties by prime superscripts. 




(M ~ ctp) dp . 


(4.1) 


As the water column in an isopycnic model is layer- 
wise discretized, the vertical integral transforms into a 
summation. In the initial, unmixed state this summation 
gives 


8 PE= S 


M k (p t 


Pd - y(P*. 


Pi) 


MSp m ~ Pd - - P«> 


(4.2) 


zations of TKE production and dissipation we have 
used. Assuming that volume is conserved during mixed 
layer entrainment also means that the surface elevation 
remains unchanged, and this allows us to evaluate the 
Montgomery potential M m of the final mixed layer as 
follows: 

M m - a„p t = M, - a,p,. (4.4) 

The PE per unit area of the final mixed layer is simply 

Pi), (4.5) 




gPE m = MJp m - pd ~ y (Pi 

and substituting from (4.3) and (4.4) above, this can be 
written 


In the final, mixed state there remains only one layer, 
having a density a m equal to the depth-weighted average 
of all the layers that have been mixed together: 

n I 

ajp m ~ Pd = 2 1 " Pd] 

k= 1 

+ [<xSP m - Pn)l ( 4 ‘ 3 ) 

Note that this expression is actually a statement of vol- 
ume conservation, which only holds for the case of a 
linear equation of state. However, the smaH error we 
make through the use of this approximation in the di- 
agnosis of the entrainment depth is of little concern, 
given the uncertainties associated with the parameteri- 


g?E m = (p m ~ pMM, - a ]Pl - 5) 


~(Pk 


- y (p m - Pn)l 


Pk) 


(4.6) 


Subtracting the expression for the potential energy of 
the initial state from this and equating the difference to 
the TKE supplied by surface stress and buoyancy ex- 
traction, 

g(P E m - PE) = gTKE, (4.7) 

leads us to an expression for the pressure at the base of 
the mixed layer following entrainment; 
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gTKE + p,(M, - a t p t - F k ) + G k - p„ 

Mn ~ ~2 

Pi + Pn ) 

Pm ~ 

(A/, - a,p, - F k ) - 

M n ~ ~(P, + Pn) 



V 


V' 




where, for convenience, we have introduced the layer- 
wise functions analogous to those of Bleck et al. ( 1 989): 


F k 


2 TT(/\ + i - Pk) and 

*= I ^ 


(4.9) 


c* = X 


Mi c (Pk+ i Pk ) 




p* 2 ) 


(4.10) 


The solution procedure for finding p m follows exactly 
that of Bleck et al. (1989). Trial values are found suc- 
cessively for each layer k = 2. . . N beneath the mixed 
layer, and the process is continued as long as the trial 
value is greater than the pressure at the upper surface 
of layer k , or the seafloor is reached. The correct value 
of p„ is the last one that meets this criterion, and once 
it has been found, updating of layer thicknesses and 
properties completes the algorithm. 


en further, with the (positive) change in PE between 
initial and final states equaling the TKE supplied over 
one time step, as before. 

While this procedure is computationally efficient, it 
suffers one potential disadvantage, in that none of the 
available PE associated w ith the unstable water column 
is dissipated, rather it all contributes to deepening of 
the mixed layer. However, in a coarse-resolution, hy- 
drostatic model any parameterization of convection is 
subject to considerable uncertainty. Given that the ul- 
timate aim of a convection routine in such a model is 
to avoid a situation that cannot be handled by the re- 
duced physics through a somewhat arbitrary reorgani- 
zation of the vertical structure, there is possibly some 
advantage in ensuring that this procedure does not rep- 
resent a sink of energy. 


c. Detrainment 


b. Convection 

In the MICOM code a separate algorithm deals with 
the case where buoyancy extraction renders the mixed 
layer statically unstable with respect to the isopycnic 
layers beneath. It is assumed that the timescale for con- 
vection is small compared to the model time step and 
thus static instability is removed by an instantaneous 
rearrangement of the appropriate model layers. Interior 
layers are incorporated successively into the mixed layer 
until the density profile is once again statically stable. 
The process is essentially one of mixed layer entrain- 
ment, but since the depth of convective mixing is di- 
agnosed in a manner that is incompatible with the theory 
outlined in the preceding section, it must be dealt with 
in a separate routine. 

We decided to exploit the similarity between the pro- 
cesses of entrainment and convection and so deal with 
them both using a single algorithm. In the case of un- 
stable stratification, deepening of the mixed layer lowers 
the PE of the water column, and the (negative) change 
in PE between mixed and unmixed states can be found 
in the manner outlined in the preceding section. The 
depth of the mixed layer following convection can also 
be found through the procedure outlined above. In the 
absence of any TKE supply, the new state will be one 
with the same PE as the initial state; that is, all the 
available PE associated with the unstable stratification 
will have been used for entrainment. With a supply of 
TKE from surface processes, the mixed layer will deep- 


In order to maintain the integrity of an isopycnic co- 
ordinate system while allowing mass to flux between 
layers one must ensure that the density of waters re- 
ceived by an isopycnic layer matches its prescribed den- 
sity. The situation of mixed layer entrainment, discussed 
in the previous sections, is not of concern in this sense 
because the mixed layer is nonisopycnic and therefore 
can receive water of any density. The reverse situation, 
that is of mixed layer detrainment, is problematic be- 
cause the density of the waters that are to be detrained 
into a particular isopycnic layer will not in general 
match the prescribed density of that layer. 

A retreat of the mixed layer occurs whenever the TKE 
generated by surface stresses is insufficient to maintain 
the current mixed layer depth in the presence of an input 
of buoyancy at the surface. In this case the TKE balance 
of Eq. (2.13) would yield a negative entrainment rate, 
but since the equation is only valid for zero or positive 
entrainment, it is used instead to diagnose a new equi- 
librium depth for the mixed layer. This is the depth for 
which the supply of TKE is exactly balanced by buoy- 
ancy and frictional dissipation, and the entrainment rate 
is therefore zero. As the mixed layer shoals to this depth 
and attains a lighter density it leaves beneath it a slab 
of water of arbitrary density, which we shall refer to as 
the fossil layer. A schematic of the relationship between 
the various densities and combinations of surface buoy- 
ancy fluxes for a detrainment scenario is presented in 

Fi g- ^ . 

A variety of solutions to the problem of how to adjoin 
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Fig. 3. Schematic of the detrainment of mixed-layer waters into an isopycnic layer as forced 
by a downward-directed surface buoyancy flux (i.e., buoyancy gain). The leftmost panel shows 
the initial state of the system with the mixed-layer density, temperature, salinity, and thickness 
denoted <r,, 0,, S { , and h ly respectively. The complementary properties of the isopycnic layer 
below the mixed layer are indicated by the k subscript. The presence of a downward surface 
buoyancy flux B m that is strong enough to suppress the TKE supplied by surface stress u] leads 
to the creation of a shallow, intermediate state layer referred to as the Caspar layer over whose 
depth scale (i.e., the Monin-Obukov length) are distributed the downward surface buoyancy 
fluxes, i.e., h g < h { . The relatively shallow Caspar layer properties are indicated in the middle 
panel by the subscript g. The difference between the shallow, intermediate state Gaspar layer 
and the deep, initial state mixed layer leads to the existence of a fossil layer in the intermediate 
state. This fossil layer is split into two sublayers: the upper fossil layer denoted with subscript 
uf and a lower fossil layer of subscript If. The final state of the system is shown in the rightmost 
panel in which the upper fossil layer adjoins with the Gaspar layer to produce the final mixed 
layer; the lower fossil layer adjoins with the isopycnic layer to produce the final state of the 
isopycnic layer. The final mixed layer properties are denoted with m subscripts and the final 
isopycnic layer properties by prime superscripts, 


the fossil layer waters to those of the isopycnic interior 
layers have been proposed, each with its own distinct 
advantages and disadvantages. One approach is to place 
all of the fossil water of arbitrary density into the iso- 
pycnic layer that most closely matches it in density (Ob- 
erhuber 1993). This has the advantage that the mixed 
layer is fully detrained and has temperature and salinity 
properties that are in accord with the buoyancy forcing 
supplied. It has the disadvantage of causing a coordinate 
drift within the isopycnic layers as they are now re- 
ceiving water masses that do not exactly match their 
prescribed densities. Another approach is to introduce 
a buffer layer sandwiched between the mixed layer and 
the isopycnic interior, which can hold the fossil layer 
waters after a detrainment event (Murtugudde et al. 
1995) thus preventing coordinate drift of the interior 
layers. This approach has the disadvantage of curtailing 
the ventilation of the isopycnic ocean interior while the 
fossil waters are held in the buffer layer. A third ap- 
proach is to split the fossil layer into a denser lower 
part, having a density precisely matching that of an 


isopycnic interior layer, and a complementary lighter 
upper part, which rejoins the mixed layer (Bleck et al. 
1992). This procedure has the advantage of ventilating 
the isopycnic interior layers with waters of the correct 
density, but the disadvantage of having to artificially 
split the heat and salt content of the fossil layer. While 
all of these approaches have limitations, these are man- 
ifestations of numerical truncation errors, which all di- 
minish in importance as the model resolution is im- 
proved. 

For our ice shelf-ocean modeling problem, the sim- 
plest approach was to follow the MICOM procedure of 
splitting the fossil layer, but we found that some mod- 
ification of the algorithm was still necessary. That used 
in the original MICOM code (Bleck et al. 1992; Sun 
1997) has some shortcomings in that under certain con- 
ditions of surface forcing it can either create an artificial 
temperature extremum or prevent detrainment proceed- 
ing at all. The problem arises because there are an in- 
finite number of ways in which the fossil layer can be 
unmixed, and selecting the most appropriate split for 
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Fio. 4. Temperature-salinity diagram under various scenarios of surface buoyancy-forced mixed- 
layer detrainment. The density of the Gaspar layer is denoted (T g and marked by t e triang e 
symbol, that of the fossil layer is denoted and marked by the circle, and that of the isopycnic 
layer is denoted <j k and marked by the square. Depending on whether heating or freshening 
dominates the surface buoyancy forcing, and whether the isopycnic layer is saltier of fresher than 
the fossil layer, gives rise to four possible density profiles. These are labeled as the , , 

■7 ” and profiles as suggested by the shape of their vertical structure and are presented in 

(a), (b), (c), and (d), respectively. 


any situation is not straightforward. While the proce- 
dures used in MICOM seem to work very well in the 
warm water sphere, for which they were originally de- 
signed, in our high-latitude polar domains their perfor- 
mance appears to degrade. 

If there were an infinite number of isopycnic layers, 
the detrainment algorithm would be straightforward. We 
would first diagnose the new thickness of the mixed 
layer from (2.13), then calculate its new temperature 
and salinity by applying the surface fluxes to this new 
layer over one time step. The fossil layer would be 
subducted directly into the interior, having properties 
identical to those of the mixed layer prior to receiving 
buoyancy input. Our detrainment algorithm mimics this 
behavior as closely as possible, given the finite number 
of isopycnic layers. 

1) MICOM ALGORITHM 

The original MICOM algorithm for detraining the 
fossil layer waters into an isopycnic layer involved sole- 
ly a temperature-based split of the fossil layer into an 
upper, warmer layer and a cooler, lower layer (Bleck et 
al. 1992). The salinity of the upper and lower fossil 
layers was not split but kept the same as the original 
fossil layer. The upper fossil layer temperature was as- 
signed precisely that of the Gaspar layer. Using these 


constraints, coupled with an equation of state, and de- 
manding that the lower fossil layer density exactly 
match that of the isopycnic layer into which the lower 
fossil layer waters are received, leads to a polynomial 
equation in one unknown, the thickness of the lower 
fossil layer. 

It was later recognized that splitting the fossil layer 
solely on temperature was problematic for polar oceans 
where the density is largely controlled by salinity var- 
iations and not temperature. Accordingly, the splitting 
of the fossil layer was adapted to involve splitting on 
both temperature and salinity (Sun 1997). In the instance 
of a density profile of shape as in Fig. 4a, the approach 
was to assign the lower fossil layer the temperature and 
salinity properties of the isopycnic layer. The upper fos- 
sil layer temperature was still assigned the Gaspar value 
as in the earlier algorithm while the upper layer salinity 
was now an unknown to be determined. As before, the 
thickness of the lower fossil layer was solved via a 
polynomial equation, which once obtained diagnosti- 
cally led to knowing the salinity of the upper fossil layer. 

In the situation that the diagnosed upper fossil layer 
salinity is less than the Gaspar or greater than the orig- 
inal fossil, as would occur in the situations depicted in 
Figs. 4b-d, then the upper fossil layer salinity is reas- 
signed to whichever of these constraining values it lies 
closer to. A new determination of the lower fossil layer 
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F[G. 5. Schematic based on a temperature-salinity diagram illus- 
trating behavior of present MICOM detrainment scheme for condi- 
tions usual to the ice shelf-ocean interaction problem. The thin dashed 
lines represent the isopycnals of the Caspar layer a*, the fossil layer 
<jj, and the isopycnic layer a k (which is the recipient of the detrained 
water). The triangle represents the Gaspar layer, which has properties 
of $ K and S g . The circle positioned with temperature 0 / and salinity 
S f represent the fossil layer properties. A built-in constraint of the 
present MICOM algorithm is that the upper fossil layer temperature 
is exactly that of the Gaspar layer. The thick, solid horizontal line of 
fixed temperature 0 U( describes this constraint. Since the unmixing of 
the fossil layer must occur along a straight line, the only possible 
detrainment lines must lie within the conic structure indicated by the 
shading. This implies the possibility of creating extremely warm tem- 
perature in the lower fossil layer or alternatively not detraining at 
all, neither behavior being desirable. 

thickness is carried out, but with the lower fossil layer 
temperature and salinity values no longer constrained 
to be exactly the same as that of the isopycnic layer. 
While the temperature and salinity of the lower fossil 
layer now differ from the isopycnic layer, the density 
of the lower fossil is still constrained to be that of the 
isopycnic layer. With these new constraints in place, the 
thickness of the lower fossil layer is once more solved 
for and the present MICOM algorithm terminates. 

In applying this algorithm to the problem of ice shelf- 
ocean interaction, situations were found where the 
above algorithm failed to detrain sufficient waters and 
the waters that were detrained were found to be unac- 
ceptably warmed. As the mixed layer in the sub-ice 
shelf cavity does not undergo a seasonal cycle, which 
would cycle the waters everywhere in the cavity through 
an entrainment and detrainment phase, regions of per- 
sistent detrainment can lead to the creation of unreal- 
istically warm water masses. Consider for instance, the 
scenario of detrainment in a polar ocean illustrated in 
Fig. 5. The fact that the present MICOM algorithm 
makes the temperature of the upper fossil layer match 
that of the Gaspar layer results in an excessively warm 
temperature for the lower fossil layer. A more natural 
split of the fossil layer would be along the direction 
“orthogonal” to the isopycnals as suggested in Fig. 5. 

Such a split avoids the generation of extreme tem- 
peratures but can only be achieved by giving up the 
constraint that the upper fossil layer temperature must 


match that of the Gaspar layer. Imposing such a con- 
straint is reasonable in a tropical ocean where density 
is dominated by temperature. In a polar ocean, where 
density is dominated by salinity, a more reasonable con- 
straint is that the upper fossil layer adhere to the salinity 
of the Gaspar layer and not its temperature. What is 
needed then is a flexible constraint on the upper fossil 
layer such that it adheres to the Gaspar temperature in 
tropical oceans but to the Gaspar salinity in polar oceans 
with a smooth transition between the two in temperate 
oceans. Such a flexible constraint is possible by intro- 
ducing the concept of an orthogonal density split of the 
fossil layer. 

2) Modified algorithm 

The main difficulty encountered with the present MI- 
COM algorithm stems form the fact that it treats tem- 
perature on a special footing as compared to salinity. 
We argue that it may be better to base a detrainment 
algorithm on density, with temperature and salinity each 
proportionally contributing to the algorithm based on 
their respective thermal expansion and haline contrac- 
tion coefficients as relevant weighting. 

We remove the above constraint that the upper fossil 
layer be assigned the Gaspar layer temperature. Instead, 
we argue that the upper and lower fossil layers can each 
take on arbitrary temperature and salinity values, pro- 
vided they are constrained by the overall minimum and 
maximum temperature and salinities present. These ex- 
trema come from the properties of the isopycnic layer, 
the fossil layer (which has the properties of the old 
mixed layer), and the Gaspar layer. The present MICOM 
scheme does not “bound” the upper fossil layer by the 
new Gaspar temperature, but rather sets it equal to it. 
We argue that the incoming buoyancy fluxes (both of 
heat and salt) are distributed over the shallower Gaspar 
layer and not the deeper, initial mixed layer. In this 
manner the fossil layer retains the initial mixed layer 
properties while the Gaspar layer has its temperature 
and salinity set by the incoming heat and salinity fluxes 
being distributed over the shorter Monin-Obukov length 
rather than the initial mixed layer depth. The new mixed 
layer temperature and salinity achieved are then used 
in defining the properties of the upper fossil layer. 

Removing the constraint that the upper fossil tem- 
perature must coincide with that of the Gaspar layer 
reintroduces an extra degree of freedom into the algo- 
rithm with the result being that there are an infinite set 
of choices for splitting the upper and lower fossil layers. 
We decide then to adhere to two guiding constraints: (i) 
the upper and lower fossil temperatures and salinities 
do not go outside the temperature and salinity extremum 
that we establish, and (ii) the lower fossil temperature 
and salinity produce a lower fossil density that exactly 
matches that of the isopycnic layer. While these con- 
straints are superficially similar to that of the present 
MICOM algorithm, the main difference lies in the 
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Taking the gradient and normalizing by density p we 
arrive at the expression 



P 


where the i and j unit vectors are aligned with the sa- 
linity and temperature axes, respectively. In the above 
equation we have made use of the conventional defi- 
nitions of the thermal expansion and haline contraction 
coefficients of seawater (albeit dimensionless for present 
purposes): 


1 da 1 da 

p d6 ^ p dS 


(4-12) 


Fig. 6. A schematic of the orthogonal density splitting of the fossil 
layer illustrated on a temperature-salinity diagram. The thin dashed 
lines are reference isopycnals showing the density of the Caspar layer 
o- g , the fossil layer <r /t and the isopycnic layer a k . The thick dotted 
line is aligned with the direction of the gradient vector (shown as the 
thick, short arrow) of the density contours and passes through the 
temperature-salinity point of the fossil layer, denoted by the circle 
symbol. The waters of the fossil layer are ultimately split into an 
upper fossil layer, denoted by the triangle, which adjoins the Gaspar 
layer to form the final mixed layer, and a lower fossil layer, denoted 
by the square, which adjoins the isopycnic layer to form the final 
isopycnic layer. The determination of the exact orientation of the 
normal vector and the amount of fossil water going into each of the 
upper and lower fossil layers is discussed in the text. 

choice of the temperature extremum. The MICOM al- 
gorithm set the upper fossil temperature to that of the 
Gaspar layer and places no constraint on the lower fossil 
temperature. Our approach is different in that we do not 
impose the upper fossil temperature, rather we solve for 
the upper and lower fossil temperatures subject to the 
constraint that no new temperature extremum are cre- 
ated. 

3) Lower fossil properties 

The new constraint we introduce to replace the pre- 
sent MICOM constraint of assigning the upper fossil 
layer temperature is to stipulate that the fossil layer is 

split based on an orthogonal density criterion. This con- 

straint says that the unmixing line for temperature and 

salinity in the fossil layer occurs along a line that is 
orthogonal to the isopycnal of the fossil layer (see Fig. 
6). Such an orthogonal split has the advantageous prop- 
erty of giving greater weight to the temperature field in 

tropical oceans while in polar oceans the salinity field 
controls the splitting of the fossil layer. 

To define the direction of such an unmixing line we 
consider the isopycnals of the temperature-salinity di- 
agram as contours of a scalar field a. For the few re- 

maining derivations in this section we treat both tem- 
perature and salinity as being dimensionless, and we do 

this so that we can compute a “gradient” in tempera- 

ture-salinity space with temperature and salinity treated 

as the independent variables of the density function. 


Given that we know the density of the fossil layer 
(which is just that of the initial mixed layer), the tem- 
perature and salinity of the fossil layer (which is just 
the temperature and salinity of the initial mixed layer), 
and the density of the isopycnic layer below, we can 
determine the orthogonal temperature and salinity of the 
lower fossil layer so that the layer achieves the same 
density as the isopycnic layer below as 

(0 lf , S v ) = - & + <*/> 3 * S * * * * />- < 4 * * - 1 3 > 

The fossil layer, prior to splitting into upper and lower 
fossil layers, has a thickness h f that is the difference 
between the initial mixed layer thickness and the shal- 
lower Gaspar layer thickness: 

h f = h { - h g . (4.14) 

The temperature and salinity of the fossil layer prior to 
its split is the initial mixed layer temperature 0, and 
salinity S { . There are six unknowns to be determined in 
the problem of splitting the fossil layer into an upper 
and lower fossil layer, namely the temperature, salinity, 
and thickness of the upper and lower fossil layers. The 
orthogonal split has set the temperature and salinity of 
the lower fossil layer reducing the unknowns to four. 
The conservation of mass for the splitting of the fossil 
layer sets the lower fossil thickness: 

h lf = h f - h u/ . (4.15) 

This leaves us with three unknowns, namely, the tem- 
perature, salinity, and thickness of the upper fossil layer. 

4) Upper fossil properties 

Two constraints we invoke for the fossil layer are that 
of overall heat and salt conservation: 

hfOf = hiAf + KAf h f s f = h lf s if + h uf S uf , (4.16) 

which in effect determines the upper fossil layer tem- 
perature and salinity. This situation now leaves us but 
one unconstrained variable, the thickness of the upper 
fossil layer. To address this last unknown we look at 
the properties of the final mixed layer, which results 
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Fig. 7. Bounded constraints on the final mixed layer temperature and salinity under various 
scenarios of heating/cooling and/or frcshening/salimfying. In all panels the Gaspar layer prop- 
erties are indicated by the triangle, the fossil layer properties by the circle, and the isopycnic 
density properties by the square. Under (a) surface heating the Gaspar temperature provides the 
upper bound on the final mixed layer temperature, while under (b) cooling the fossil layer 
temperature provides the upper bound, (c) Under conditions of surface freshening the Gaspar 
layer salinity provides the lower bound on salinity, while (d) under salinifying conditions, the 
fossil layer provides the lower bound on salinity. 


when we adjoin the Gaspar layer with the upper fossil 
layer. Invoking heat and salt conservation for this mix- 
ing process produces the constraints 

hj m = h g B x + h uj 6 u{ h m S m = h g S g + h uf S u/ , (4.17) 

where from mass conservation the thickness of the new 
mixed layer is simply 

h m = K + h„. (4.18) 

It appears that we have gained little as we have added 
two more equations with two more unknowns, namely, 
the final temperature 0 m and salinity S m of the final 
mixed layer. As we shall see, however, we can in fact 
ultimately combine these equations to solve for the 
thickness of the upper fossil layer. This amounts to mak- 
ing decisions about the final values that the mixed layer 
temperature 6 m and salinity S m will achieve under var- 
ious scenarios of either net heat gain or loss and sim- 
ilarly for freshwater gain or loss with the objective of 
maximizing the amount of detrainment. Specifying both 
the new mixed layer temperature and salinity now puts 
us in the situation of having one too many equations in 
relation to the number of unknowns. This situation 
means that we will obtain two physically valid solutions 
for the upper fossil layer thickness and we will take the 
thicker of these as that will represent the most con- 
strained or conservative solution. We now show' these 
two solutions. 

Combining the heat conservation equation (4.16) for 


the fossil mixed layer with the heat conservation equa- 
tion (4.17) for the final mixed layer leads to an ex- 
pression for the thickness of the upper fossil layer: 


h f (0j - 0 (f ) + - 9J 

(0 m - e lf ) 


(4.19) 


Similarly, combining the salt conservation equations we 
obtain 


h f (S f ~ S lf ) + h 8 (S H - SJ 
(S m - V 


(4.20) 


We cannot evaluate these expressions until we come to 
a decision as to the final mixed layer temperature and 
salinity. Our objective is to make the upper fossil layer 
as thin as possible so as to detrain as much water into 
the isopycnic layer via the lower fossil layer as is pos- 
sible, without violating the constraints we impose. This 
objective will be met by maximizing the difference be- 
tween the orthogonal temperature and salinity and the 
final mixed layer temperature and salinity, respectively, 
while adhering to the overall temperature and salinity 
constraints that ensure we do not create new temperature 
or salinity extrema. It turns out that this stipulation of 
maximum difference depends upon, for temperature, 
whether we are cooling the Gaspar layer or heating it. 
Analogously for salinity, we need to know whether the 
Gaspar layer is freshened or salinified. These situations, 
depicted in Fig. 7, give us a physical interpretation that 
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guides the assignment of the final mixed layer temper- 
ature and salinity. Briefly stated, if the Caspar layer is 
heated, then the final mixed layer temperature will at- 
tempt to be the Gaspar layer temperature (Fig. 7a), while 
if the Gaspar layer is cooled, the final mixed layer tem- 
perature will be the warmer fossil layer temperature 
(Fig. 7b). Analogously for salinity, under conditions of 
freshwater input to the Gaspar layer via buoyancy forc- 
ing then the final mixed layer salinity will attempt to 
lower to that Gaspar value (Fig. 7c), and vice versa for 
salt input to the Gaspar layer the final mixed layer sa- 
linity will be the fresher fossil layer salinity (Fig. 7d). 
In this manner the maximum amount of detrainment, 
while not creating new temperature or salinity extre- 
mum, is achieved. 

The four scenarios outlined above give four different 
physically meaningful solutions to the upper fossil layer 
thickness solutions (4.19) and (4.20). Under the case of 
heat input (Fig. 7a) at the ocean surface the upper layer 
fossil thickness is 


0 f - 0 lf 


(4.21) 


which is a well-behaved solution as we can formally 
guarantee the desired bounded property of 0 ^ h uf ^ 
hj. Under the case of heat extraction (Fig. 7b) the so- 
lution is 


0 f - 0 g 

h ■' - *' - 


(4.22) 


which becomes a well-behaved solution as we restrict 
the permissible solutions to the desired range of 0 ^ 
h uf < h f . Under the case of freshwater input (Fig. 7c) 
the solution is analogous to (4.21) with salinity substi- 
tuted for temperature. Under the case of salt input (Fig. 
7d) the solution is analogous to (4.22) again with salinity 
substituted for temperature. 

Of the two temperature-based upper fossil layer thick- 
ness solutions (4.2 1 ) and (4.22), only one can be realized 
at any given instant in time as the Gaspar layer is either 
being heated or cooled but not both simultaneously. 
Similarly, only one salinity-based thickness solution is 
physically relevant at a given instant. This finally leaves 
us with one temperature-based upper fossil layer thick- 
ness estimate and one salinity-based one. Again, we 
choose the larger upper-fossil thickness as the one to 
represent the most constrained solution, whether guided 
by temperature or salinity extremum considerations. 
This then is the solution of the upper fossil layer thick- 
ness. All other dependent quantities of the upper and 
lower fossil layers are now determined and the algorithm 
is complete. 

The advantage in this approach is clear in that we can 
always detrain while not altering the isopycnic density 
of our interior layers, and that we treat tropical water 
masses and polar water masses on an equal footing with 
a smooth transition between the two types because our 


split is on Gaspar density and not on Gaspar tempera- 
ture. 

d. Surface buoyancy fluxes 

Because the base of the ice shelf is a thermodynam- 
ically active boundary, the waters directly in contact 
with it will be in a continual state of adjustment toward 
conditions of thermodynamic equilibrium. We use this 
fact to diagnose the surface heat and freshwater fluxes 
that are applied to the mixed layer beneath the ice shelf. 
The procedure is described in detail by Holland and 
Jenkins (1999). We use their three-equation formulation, 
which explicitly accounts for the differing diffusivities 
of heat and salt in the inner part of the ice-ocean bound- 
ary layer, the molecular sublayer. Heat conduction into 
the ice shelf is parameterized in a way that accounts 
approximately for both vertical diffusion and constant 
vertical advection within the ice. The melting and freez- 
ing that occur at the ice base are the physical processes 
responsible for this vertical advection. 

Throughout the model domain, whenever we compute 
the response of the mixed layer to surface forcing, we 
always convert from potential to in situ temperature 
before undertaking the calculation. The altered in situ 
temperature is then converted back to potential tem- 
perature for use in all other model calculations, for ex- 
ample, pressure gradients or mixed layer entrainment/ 
detrainment. This back and forth transformation be- 
tween potential and in situ temperature, which is carried 
out in double-precision arithmetic, ensures the correct 
application of surface heat fluxes, regardless of the sur- 
face pressure, and does not result in any net drift in 
temperature. 

e. Thermodynamic initialization 

The waters throughout the cavity are first assigned 
properties through lateral extrapolation of values found 
at the ice front. If the assigned properties of the mixed 
layer are far removed from a condition of thermody- 
namic equilibrium with the ice, there will be an initial 
spike of high melting, which could potentially cause a 
numerical instability. Therefore, before taking the first 
model time step, we relax the temperature and salinity 
of the mixed layer beneath the ice shelf to a state of 
thermodynamic equilibrium. 

We first calculate the in situ freezing point 7} based 
on the initial assigned temperature T a and salinity S a 
and the pressure at the ice shelf base p x . 

T f = + A 2 + A 3 / 7 , , (4.23) 

where A, 3 are empirical constants (Millero 1978). The 
adjustments we make mimic the effects of melting or 
freezing on the mixed layer properties, so the temper- 
ature and salinity are constrained to evolve along a 
straight line in T/S space (Gade 1979) given approxi- 
mately by 
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dT L f 

ds ~ c F s a * 


(4.24) 


where L f is the latent heat of fusion and c p is the specific 
heat capacity of seawater at constant pressure. The sa- 
linity of the mixed layer waters, once they have reached 
the freezing point, can then be expressed in terms of 
the initially assigned conditions: 





(4.25) 


As a final adjustment, the mixed layer temperature is 
recomputed with the salinity from (4.25) using (4.23). 
We emphasize again that the above expressions are all 
in terms of the in situ temperature, T t rather than the 
potential temperature, 0. 


5. Interactions with topography 

a . Layer discontinuity 

As a result of the horizontal discretization, surface 
and bottom topography acquires a stepped structure 
within MICOM. This raises the possibility of a layer 
becoming discontinuous, if it is thinner than a particular 
topographic step is high. In our sub-ice shelf applica- 
tion, this is a particularly strong possibility, because 
both sea surface and seabed have topography. In par- 
ticular, the vertical front of the ice shelf introduces a 
surface step of up to several hundred meters, making it 
likely that the mixed layer is one of the layers containing 
a discontinuity. A particular example of this was earlier 
illustrated in Fig. 1 in the case of the mixed layer. 


1) Advection 

In principle, there is no difficulty in estimating the 
mass flow that should pass between grid points that are 
separated by a topographic discontinuity. The calcula- 
tion of the Montgomery potential gradient is unaffected 
by such a discontinuity and will correctly drive flow in 
the direction that is energetically favorable, that is, 
dense layers will be forced down bottom slopes, light 
layers up surface slopes. However, this does not in itself 
ensure that mass and tracers can be advected across a 
discontinuity. 

The cause of the difficulty arises from the fact that 
MICOM is discretized on an Arakaw'a C grid (Arakawa 
and Lamb 1977) in which the thickness and other scalar 
properties of layers are defined at grid points horizon- 
tally staggered with respect to the points at which the 
velocity components are defined. A similar problem 
would also arise if the discretization was performed on 
the Arakawa B grid. In order to compute the flux of 
mass or tracer, using (2.5) or (2.6), from one scalar point 
to a neighboring one, a layer thickness must be assigned 
to the velocity point that lies between the tw'o neigh- 
boring scalar points. The manner in which this assign- 


ment is made in the original MICOM code and in our 
modified code is illustrated in Fig. 8. 

Consider the specific example of a two-layer fluid 
overlying a continental shelf and adjacent deep ocean. 
The surface of the top layer has a pressure denoted p x 
and the bottom of the bottom layer has a pressure /> 3 
while the pressure along the mutual interface between 
the two layers is denoted p 2 . We horizontally discretize 
this fluid environment by considering two “scalar” col- 
umns, labeled column a and c, with column a being 
representative of the shallow layer of fluid on the con- 
tinental shelf and c representing the deep-ocean column. 
A scalar column here refers to a column that contains 
vertically coaligned scalar variables of the Arakawa C 
grid; a “vector” column refers to vertically coaligned 
velocity points. Vector columns occur halfway between 
neighboring scalar columns. In Fig. 8, halfway between 
the columns a and c we have such a vector column, 
labeled b. The assignment of layer thicknesses at the 
vector column requires that we establish the pressure p 2 
at the position of the vector column. Referring to Fig. 
8, we use X to mark where this pressure is assigned to 
occur for the original MICOM code and our modified 
code. In the original MICOM code, p 2 coincides exactly 
with Pi at a vector column. The thickness of the bottom 
layer at the vector column is then zero; that is, the thick- 
ness is evaluated as p 3 — p 2 — 0. Tn the modified code, 
p 2 is assigned to occur higher up in the vector column 
b such that p 3 - p 2 ^ 0, and thus a nonzero thickness 
of the bottom layer at the vector column is prescribed. 

To understand better how this problem can at all oc- 
cur, it is first necessary to realize that the vector column 
b has a total column thickness defined by the difference 
of the deeper of the neighboring surface elevations at 
columns a and c and the shallower of the two seabed 
elevations at a and c. Given such a thickness of the 
vector column b there still remains a degree of freedom 
in assigning the thickness of the upper and lower layer 
within this vector column. It is that assignment that is 
the key difference between the original MICOM code 
and our modified code. 

The algorithm usecTTn the original MICOM code as- 
signs zero layer thickness to the bottom layer at the 
vector column b. This is so because it uses the maximum 
“seafloor depth” to evaluate the layer thickness at vec- 
tor points. Visually one can see this by inspecting the 
low'er-Ieft panel of Fig. 8 and noting that the original 
MICOM assigns the position of X, which represents the 
value of p 2 at vector column b, as the average height 
between p 2 at column a and the p 2 at column c. Mass 
flux in the bottom layer from column a to column c 
would then be forced to zero because the bottom layer 
has zero thickness at the vector column. This would 
lead to the entrapment of dense water on the continental 
shelf. Tn the ice shelf-ocean context, the upside-down 
version of this problem w ould occur and buoyant water 
w'ould be trapped beneath the deeper parts of the ice 
shelf base with no possibility of floating upward to lower 
elevations. 
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Fig. 8. Illustration of the problem of a dense lens of fluid over a continental shelf being 
artificially confined to rest on the shelf while in fact the dense plume should flow off the shelf 
into the deeper waters. The top, centered panel is intended to show a continuum representation 
of a situation of a two-layer fluid overlying a continental shelf and slope. In the present MICOM 
algorithm, shown in the leftmost lower panel, the discretization of this scenario leads to a situation 
whereby there is no "connection" between the fluid overlying the shelf and that in the deeper 
ocean; consequently there is no flow off the shelf. In the modified algorithm, shown in the rightmost 
lower panel, there is a connection between shelf and slope waters, and thus the dense lens of 
water on the shelf can flow off the shelf and into the deep. The key difference is that the modified 
algorithm uses the concept of a shelf horizon depth instead of a seafloor depth to compute the 
thickness of the layers at velocity points on the Arakawa C-grid stencil. See the text for further 
details including an explanation for the various symbols used in the figure. 


We introduce the term shelf horizon depth to mean 
the shallower of the bottom ele vation of the two columns 
a and c, which in this instance occurs at column a. We 
assign layer thicknesses at vector column b by using the 
shelf horizon depth concept to produce a nonzero thick- 
ness for the bottommost layer, in a physical situation 
where that is an appropriate assignment. This can be 
seen visually in the lower-right panel of Fig. 8 where 
the value of p 2 and in column c is not allowed to 
fall below the shelf horizon depth. Our algorithm works 
to assign the positiort of the X higher up in the vector 
column b than in the case of the original MICOM code. 
Our approach is, however, similar to the original MI- 
COM scheme in that is also ensures that the sum of all 
the layer thicknesses does not exceed the total water 
column thickness at the vector column. The distinguish- 
ing feature is that any layer with nonzero thickness in 
either of the neighboring scalar columns a or c is now 
assigned a nonzero thickness at the vector column b. 
The result is that gravity currents flow freely both for 
the case of a dense plume descending a continental slope 
and that of a buoyant plume ascending the base of an 
ice shelf. 


2) Diffusion 

While advection is allowed to proceed across a layer 
discontinuity, diffusive fluxes are modified according to 
a horizontal “line-of-sight” criterion. If a given layer 
abuts other isopycnic layers, the diffusivities and vis- 
cosity are simply set to zero, whereas if it abuts a solid 
wall (either ice shelf or seabed), scalar diffusivities are 
set to zero, while the standard viscosity is used. In this 
latter case velocity boundary conditions are applied con- 
sistent with either free slip or no slip at the solid wall. 

If a layer is not completely disjointed by topography, 
but a partial overlap exists between neighboring scalar 
grid points, normal diffusive exchange is permitted but 
fractionally weighted according to the proportion of the 
layer that is in contact with the solid wall. In this manner 
there is a smooth transition of the diffusive fluxes if the 
discontinuity in a layer should gradually appear or dis- 
appear. 

b. Layer thickness diffusion 

As the horizontal resolution of present-day ocean gen- 
eral circulation models is often inadequate to resolve 
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the detailed structure of baroclinic eddies, a parameter- 
ization of their impact is utilized. In the isopycnic frame- 
work this is layer thickness diffusion, and is represented 
by the term on the right-hand side of (2.5). This term 
acts to flatten our isopycnals, so as to make them co- 
incide with geopotential surfaces, and mimics the way 
mesoscale eddies remove potential energy from the 
ocean interior. In the present MICOM code, the flatten- 
ing of the isopycnals is prohibited once this process 
would force layer interfaces to go deeper than the bed- 
rock. In the situation of an imposed surface topography, 
the algorithm is modified so as to also prevent isopyc- 
nals passing upward shallower than the base of the ice 
shelf. 

6. Interfacial friction 

The horizontal momentum equations (2.8) include a 
stress term rjf that describes the vertical flux of hori- 
zontal momentum between neighboring layers. These 
viscous stresses can occur on the interface between any 
two adjacent layers or between a layer and an adjacent 
solid boundary. Such interfacial friction is commonly 
parameterized in ocean models using a quadratic drag 
law of the form 

Tf = p„c d |Av;-|Avi, (6.1) 

where c d is a nondimensional and empirical drag co- 
efficient, typically of order 10 3 . The velocity jump 
across the lower 4- and upper - interfaces bounding 
layer k are defined 

Av* = v*+, - v k and Av^ = v k - v k _ r (6.2) 

In MICOM currently only bottom friction is explicitly 
included. The velocity at the seabed is assumed to be 
zero, so the velocity jump appearing in the drag law 
(6.1) is simply the layer velocity. For simplicity, the 
drag law is linearized through the use of an effective 
drag coefficient c d , based on the velocities from the 
previous time step, denoted by the “old” superscript: 

c\i ~ p„c d \ Av; | old . (6.3) 

We note that this linearization renders the computation 
of interfacial stress susceptible to numerical instability 
under certain flow conditions, but we do not offer an 
alternative. The next step in the MICOM algorithm is 
to define a viscous thickness scale (A p) v , typically 
equivalent to 10 m of water, over which the frictional 
stresses are actually applied. Any fluid, of whatever den- 
sity, found within this distance of the seafloor is sub- 
jected to drag. The viscous stress applied to each iso- 
pycnic layer is linearly weighted by the fraction of that 
layer residing within the bottom (A/?)„ of the water col- 
umn. Although the velocity vector of any layer that is 
subject to bottom drag will turn away from that of a 
purely geostrophic flow, the above formulation cannot 
provide a general parameterization of the bottom Ekman 


layer, (Ekman 1905), even if there is sufficient model 
resolution near the seabed. 

With the addition of static sea surface topography, 
surface friction could be introduced in a manner anal- 
ogous to that described above for bottom friction. The 
algorithm would be more straightforward, because the 
mixed layer, aside from being the only nonisopycnic 
layer, is distinct from the other layers in that it is the 
only one that has a minimum specified thickness (typ- 
ically equivalent to 10 m of water). As a consequence, 
only the mixed layer would experience the surface stress 
generated at the ice shelf base. For the sake of com- 
pleteness, we mention that in the open ocean, away from 
the ice shelf base, the mixed layer experiences a surface 
stress due to the presence of sea ice and the atmospheric 
winds. 

Aside from surface and bottom stresses, we also wish 
to consider lay er-to- layer interfacial stresses in the 
ocean interior that arise both from small-scale mixing 
processes and as” a result of the breaking of the internal 
wave field. At the interface between any two adjacent 
layers we impose the quadratic drag law, but with a drag 
coefficient that is much reduced compared with that used 
for the seabed or surface friction. Our general formu- 
lation of the viscous stress is then Eq. (6.1) applied to 
every layer interface, but with the constant interface 
drag coefficient, c d , replaced by a variable drag coef- 
ficient, c" d : 

c'J = c d [[e~* Pi ~ Px l /(Ap)v + e~^ Pk ~ pb ^ ap)v ] + 0.001}. (6.4) 

This formula has the desired effect of making the ef- 
fective drag coefficient, c d , equal to the standard value, 
c rf , near the ocean bottom, p = p h < and surface, p = p ,, 
and allowing it to decay gradually to a value two orders 
of magnitude smaller far away from the boundaries. 

For an TV-layer ocean model, there are N horizontal 
velocities. Our algorithm for the application of inter- 
facial stresses is to simultaneously solve (6.1) for all 
the layer velocities. Applying the boundary conditions 
of zero velocity at the seabed and ice shelf base, w F e 
arrive at a system of N nonlinear equations. We linearize 
the drag law by adopting an effective drag coefficient 
based on velocities at the previous time step, in analogy 
with (6.3). As the net viscous stress acting on any given 
layer is a function of only the velocities of that layer 
and of the layers directly above and below it, we then 
arrive at a system of N linear equations in tridiagonal 
form. These are conveniently solved using a standard 
implicit technique, such as Gaussian elimination. It is 
also noted that this formulation is capable of producing 
an “Ekman spiral,” given sufficient resolution within 
the physical boundary layer. 

7. Model application 

a. Configuration 

We have applied the model to a 10° long X 10° lat 
X 1000-m-deep box with an idealized ice shelf floating 
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F : [(i . i ). Idealized box domain with a lloaiing ice shelf in the southern 
half of the domain (north is to the right). The iee shell front has a 
uniform thickness of 200 m while at the grounding line iU\. where 
the ice shelf meets the bottom topograph)- ) the thickness extends to 
700 m There is a uniform slope between the ice shell tront and the 
grounding line. The model domain extends south-north over 10 ot 
latitude, or equivalent!), about 1000 km: the domain extends west- 
east over 10° of longitude, or equivalently, about 2.^0 km. The max- 
imum depth of 1 0(10 m ot the water column occurs in the open-ocean 
part of the domain: the minimum depth of 200 m occurs beneath the 
ice shelf along the grounding line. 


in the southern half of the domain (sec Fig. 9). The 
central latitude, which coincides with the ice shell front, 
is 75 C S, so the horizontal dimensions of the domain are 
approximately 250 km in width east-west by 1000 km 
in length north-south. In thickness, the ice shelf itself 
is wedge shaped with the thickness increasing linearly 
from 200 m at the front to 700 m at the southern bound- 
ary of the domain and being zonally constant. The grid 
is isotropic w ith the grid cells having dimensions equiv- 
alent to 0.5° of longitude. This gives a horizontal res- 
olution that ranges from approximately 10 km in the 
southern part of the domain and nearly 20 km in the 
northern. We divide the water column north of the ice 
shelf front into 1 1 layers (10 isopycnie layers plus the 
mixed layer) of equal thickness and extend the layer 
interfaces horizontally below the ice shelf. The potential 
temperature in all the layers is initially “1.8 C, while 
the salinity increases in equal steps from 34.4 psu in 
the open ocean mixed layer to 34.8 psu in the bottom 
layer. Beneath the ice shelf the mixed layer first takes 
on the properties of the isopycnie layer it displaces, and 
is subsequently relaxed to the freezing point, as de- 
scribed in section 4e. Restoring conditions arc used for 
the temperature and salinity along the northern wall so 
as to allow the model to evolve toward a steady state. 
Without such restoring, the only thermodynamic forcing 
w 7 ould be that occurring at the ice shelf base and this 
in itself would lead to a gradual cooling of all of the 
waters in the domain. 

b. Initial tests 

The model setup described above, but with no ther- 
modynamic initialization of the sub-ice shelf mixed lay- 
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Fit;. 10. Homogeneous ocean with no ice shelf base thermods- 
namical forcing, (a) Depth-averaged How velocities (cm s * ) at the 
end of the 20-yr model simulation. The position ot the ice shelf Iron! 
is at 7S S. (b) Time series over the last 10 yr ot the model integration 
of the domain integrated zonal and meridional depth-averaged flow s. 
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er or any subsequent thermodynamic interaction be- 
tween ice and ocean, has an exact solution ot zero mo- 
tion for all time. Given that wc know of no analytical 
solutions for buoyancy-driven circulation within a cav- 
ity, reproduction of such a motionless state provides us 
w ith our one chance of verifying at least some of our 
model code against a known solution. We show here 
the results of two such tests, one in which the ocean 
was completely uniform, with a salinity of 34.8 psu in 
all layers, and one in which the ocean was stratified in 
the manner described above. 

Depth-averaged velocities at the end of a 20-yr sim- 
ulation of the homogeneous ocean are shown in Fig. 
10a. Peak values of the order 10 s cm s 1 occur along 
the first row of grid points below the ice shelf, while 
the domain average is four orders of magnitude smaller 
(Fig. 10b). There is a linear growth in the magnitude 
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Fig. II, Strati lie J ocean with no ice shelf base thermodynamical 
forcing, (a) Depth-averaged flow velocities (cm s' 1 ) at the end of 
the 20- yr model simulation, (b) Time series over the last 10 \ of the 
model integration of the domain-integrated zonal and meridional 
depth-aseraged flows. 



Fig. 12. Plan view of ice shelf basal melt rate (cm yr ' ). The 
northern half of the domain represents the open ocean and the south- 
ern half the ice shelf base. The red-shaded areas represent melting 
zones and the blue represent freezing zones. It is in the later where 
marine ice accumulates at the ice shelf base. The mean melt rate is 
0.5 cm yr 

cause discontinuities in the mixed layer the front can 
be preserved. The density gradient within the mixed 
layer then exactly balances the gradient in Montgomery 
potential, and the motion remains insignificant. As the 
horizontal grid resolution increases toward the south, 
the vertical steps in the topography become smaller. 
South of 77.2°S the surface steps are less than 10 m, 
so the mixed layer is continuous in this part of the 
domain, and the density fronts within it cannot be main- 
tained against diffusion. The exact balance of forces is 
disrupted and zonal bands of motion develop at the 
points where isopycnic layer interfaces intersect the 
mixed layer. 

The above results illustrate that genuine noise in the 
depth-averaged velocity fields is stable at the level of 
10 * cm or less. This gives us confidence that we 
have implemented the surface pressure field and asso- 
ciated model modifications in a manner that has not 
caused any spurious forcing on the ocean. 


of the zonal component with time, but for any con- 
ceivable integration time this noise will remain incon- 
sequential. 

Analogous results for the stratified case are shown in 
Figs. 1 la and 1 lb. In this case peak velocities of the 
order I mm s' ? occur along three bands in the most 
southerly part of the cavity. Once again domain-aver- 
aged values are four orders of magnitude smaller. While 
there is no growth in magnitude over the final 10 yr ot 
the integration, and the general noise level is insignif- 
icant, the same cannot be said of the peak values. How- 
ever, these velocities have a physical origin. They reflect 
the fact that w ithin the mixed layer, isopycnic surfaces 
arc no longer horizontal, hut are effectively vertical. In 
our discretized space this manifests itself as a sharp 
density front in the mixed layer wherever its base in- 
tersects an isopycnic layer interface. Discretization also 
leads to steps in the surface topography, and where these 


c. Fid! model 

When the stratified domain is forced by the ther- 
modynamic interaction between the mixed layer and the 
ice, a circulation pattern emerges that, after 20 yr, leads 
to net melting of 0.5 cm vr ’ averaged over the ice shelf 
base. There is spatial variability (Fig. I2j with melting, 
at peak rates of around 10 cm yr concentrated at the 
southern boundary and broad areas of more gentle melt- 
ing and freezing to the north. The zone of freezing oc- 
cupies the central and western portion of the ice shelf 
base with peak rates of around -0.5 cm yr 1 occurring 
in the northwestern corner. Overall, the thermodynamic 
interaction with the ice shelf supplies a surface fresh- 
water flux of 22 m 3 s', and the resulting thermohaline 
circulation drives a transport of 0.06 Sv into and out of 
the sub-ice shelf cavity. The circulation beneath the ice 
shelf is cyclonic, with inflow concentrated in the east 
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Pic*. 13. Fast-west vertical transects across the ice shelf front. The 
view is from within the ice shelf cavity and looking northward, (a) 
The meridional velocity (cm s } shows northward flow as red shad- 
ing in the upper left while the lower right shows southward flow as 
blue shading, (b) The potential temperature C showing the presence 
of supercooled fSW as blue shading near the base of the ice shelf 
and somewhat concentrated to the western side of the front. The 
deeper waters continue to hold their initially assigned temperatures 
of - F8°C and are shown as red shading. 

and outflow in the west. A vertical section of meridional 
velocity across the ice front (Fig. 13a; shows a level of 
zero motion sloping down from east to west, consistent 
with a picture of deep inflow, upwelling within the cav- 
ity and outflow along the ice shelf base. Peak velocities 
reach only a few millimeters per second. Near the ice 
shelf base the outflowing water has a temperature of 
around ~2°C (Fig. 13b) and is, thus, classed as ISW. 
Deeper in the water column, much of the outflow has 
been unaffected by the ice shelf and has retained its 
initial temperature. 

The cyclonic gyre beneath the ice shelf shows up in 
the depth-mean flow (Fig. 14a). Viewed in this way, the 
gyre is intensified to the west and is centered just south 
of the ice shelf front. Although much of the depth-mean 
flow appears to be steered along the ice front, the step 
change in water column thickness does not form a com- 
plete barrier to the circulation (Fig. 14b). There are 
broad regions at the east and west where the velocity 
vectors have a significant component across the ice shelf 



Longitude 


Fig. 14. Bumtropic characteristics of the domain and the flow field 
at the end of the 20-yr model simulation, (a) The depth-averaged 
velocity field (cm s ’) for the stratified ocean with active thermo- 
dynamical forcing at the ice shelf base, (h) The contours of back- 
ground "burotropic” potential vorticity 1 1 /( m s > X 10 ' ] emphasize 
the influence exerted by the ice shelf basal topography. 


front. This is not surprising considering that the only 
forcing on the circulation is the tilting of isopycnic sur- 
faces by the freshwater flux at the ice shelf base. The 
forcing on an individual model layer is largely inde- 
pendent of that on the other layers, and a nonzero depth- 
mean flow emerges only because of the influence of 
friction and the presence of the surface topography it- 
self. This point is well illustrated by the flow of the 
mixed layer (Fig. 1 5a), which shows no hint of the gyre, 
but instead illustrates behavior we would associate with 
a buoyant plume ascending the ice shelf base. Because 
the mixed layer remains thin over much of the cavity 
(Fig. 15b), the dominant balance of forces is between 
buoyancy and friction, and deflection by the Coriolis 
acceleration is relatively minor. The mixed layer flow 
is largely sustained by the entrainment of fluid along 
the southern boundary, where there is strong horizontal 
divergence, and the same occurs to a lesser extent along 
the eastern boundary. This pattern of entrainment gives 
rise to the regions of strong melting described above. 
As the waters of the mixed layer ascend the ice shelf 
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Fir. 15. Mixed layer characteristics at the end of the 20-yr model 
simulation, (at Mixed layer How Held (cm s ') illustrating the north- 
westward buoyam-plume-like flow regime of the mixed layer cir- 
culation beneath the ice shelf, (b) Mixed layer thickness (in) showing 
a rather thin mixed layer beneath the ice shelf that thicknens in the 
open-ocean portion of the domain. 


base the decreasing pressure leads to supercooling and 
the onset of freezing. 

(I. Discussion 

The sub-ice shelf circulation presented above is 
broadly similar to that described in earlier studies. The 
concept of a buoyant flow of ISW ascending the ice 
shelf base after being initiated by mixing and melting 
in the deepest parts of the cavity, much like the behavior 
of our mixed layer, was first introduced bv Robin ( 1 979). 
The idea was subsequently developed by MacAyeal 
(1984; 1985b) and formed the basis of the one-dimen- 
sional models of Jenkins ( 1991 ) and Jenkins and Bom- 
boseh (1995). Helmer and Olbers (1989) found similar 
behavior in a two-dimensional, vertical plane model. 
With the advent of three-dimensional modeling, atten- 
tion shifted away from the overturning circulation. De- 
termann and Gerdes (1994) and Grosteld et al. (1997) 
concluded that the circulation was dominated by a 


strong depth-averaged horizontal flow. In applications 
to idealized domains, they found that the depth-inte- 
grated transport formed a cyclonic gyre in the cavity, 
qualitatively similar to the gyre indicated by the depth- 
averaged flow' vectors in Fig. 14a. However, the results 
we have presented here are quantitatively quite different 
to those of all the earlier studies. The rates of melting 
and freezing and the strength of both the horizontal and 
overturning circulations are all an order of magnitude 
smaller in our model results. 

The likely cause of at least some of these differences 
is the differing parameterizations of vertical diffusion. 
In a model where melting at the ice-ocean interlace is 
the only source of forcing, the transport of heat to the 
ice shelf base is of paramount importance. II vertical 
heat transfer is easy, melting will be rapid and the cir- 
culation will receive strong forcing. However, there is 
a natural brake on this system in that strong melting 
suppresses vertical mixing. This brake is incorporated 
into our model through the use of the TKE budget of 
the mixed layer to diagnose entrainment rates. In con- 
trast, the three-dimensional models of Determann and 
Gerdes (1994) and Grosfeld et al. (1997) and the two- 
dimensional model of Hellmer and OTbers (1989) used 
constant eddy diffusi vities to estimate vertical fluxes, 
and the one-dimensional models of MacAyeal (1985b), 
Jenkins (1991), and Jenkins and Bdmbosch (1995) all 
used a simple entrainment scheme that did not consider 
the surface freshwater flux. The impact of the negative 
feedback of melting on vertical heat transfer is evident 
in Fig. 12. Only in the vicinity of the grounding line, 
where entrainment of warm water is driven by horizontal 
divergence in the mixed layer flow, can melting in ex- 
cess of I cm yr _ 1 be sustained. Introducing an additional 
source of energy for mixing, such as strong tides, to the 
model would not necessarily lead to greater or more 
widespread melting. While a greater supply of TKE 
should increase the overall depth of the mixed layer, 
continuous entrainment would still only occur at those 
locations where it is required to maintain the thickness 
against thinning driven by horizontal divergence of the 
mixed layer flow. 

Another difference between our model and the three- 
dimensional models of Determann and Gerdes (1994) 
and Grosfeld et al. (1997) is the manner in which heat 
and salt transfer at the ice shelf base is parameterized. 
The various formulations of this process that have been 
used in models of sub-ice shelf circulation are discussed 
by Holland and Jenkins ( 1999). That used in our model 
makes the rate of heal and salt transfer dependent on 
the friction velocity, while that used by both Determann 
and Gerdes ( 1994) and Grosfeld et al. (1997) effectively 
assumes a constant friction velocity. The two formu- 
lations would give equal melt rates only if the mixed 
layer velocity were about 35 cm s 1 (Holland and Jen- 
kins 1999). Our weaker flow thus gives much weaker 
forcing, which once again contributes to the weaker 
flow. 
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The results we have presented here differ from those 
of earlier three-dimensional modeling studies in another 
important respect. The major conclusion of both De- 
termann and Gerdes (1994) and Grosfeld et al. (1997) 
was that the circulation beneath the ice shelf remains 
independent of that in the open ocean unless “////” 
contours cross the ice shelf front. That conclusion is not 
supported by our results. Our model forcing is confined 
beneath the ice shelf, but the resulting circulation ex- 
tends over the entire domain and drives exchange be- 
tween the cavity in the south and the open ocean to the 
north. 

8. Closing remarks 

The study of ocean circulation beneath ice shelves 
has, of necessity, progressed largely without the aid of 
the direct observational evidence normally taken for 
granted in oceanography. There exists a growing da- 
tabase of ship-based observations made along ice fronts, 
but these features themselves exert a strong influence 
on the water column. They guide a vigorous coastal 
current, are the sites of wind-forced upwelling and 
downwelling (van Heijst 1987) and possibly also of 
thermohaline convection (Foldvik and Kvinge 1974), 
and represent sharp changes in water column thickness 
that induce rectified tidal flows (MacAyeal 1985a; Mak- 
inson and Nicholls 1999) and may act as a barrier to 
some of the flow in the open ocean, such as that forced 
by the wind (Grosfeld et al. 1997). Using ship-based 
observations to study sub-ice shelf processes is like try- 
ing to understand continental shelf processes based only 
on data collected along the shelf break. While gross 
budgets and water mass transformations can be inferred, 
only theory, supported by a few isolated observations 
(Jacobs et al. 1979; Nicholls and Makinson 1998), can 
guide us to the important mechanisms that effect those 
changes. Yet here we encounter another difficulty. The 
most tractable theoretical problems in oceanography are 
those in which the ocean can be assumed to be ho- 
mogeneous. Although the ocean beneath an ice shelf 
may appear to be weakly stratified, spatial variations in 
the stratification are what drives the flow. The assump- 
tion of a homogeneous water column yields the null 
result discussed in section 7b. 

Starting with the work of Determann and Gerdes 
(1994), recent efforts to understand processes beneath 
ice shelves have involved the application of suitably 
modified ocean general circulation models to the prob- 
lem. The model presented in this paper is the first to 
employ isopycnic coordinates. Such a choice of vertical 
coordinate is natural for the density-driven flows of in- 
terest in this study. It ensures optimum resolution of the 
water column and allows the introduction of entirely 
arbitrary topography while avoiding the problems en- 
countered in geopotential coordinates with up- or down- 
slope flows or the difficulties encountered in terrain- 
following coordinates with estimating pressure gradi- 


ents. Precise control over diapycnic diffusion and the 
exact conservation of potential vorticity are also more 
easily implemented in isopycnic coordinates. The one 
drawback is the requirement for a nonisopycnic layer, 
which can admit buoyancy forcing, at the top of the 
water column. Here some of the advantages outlined 
above are lost, and the use of an embedded mixed layer 
model adds complexity, particularly when dealing with 
the surface topography imposed by an ice shelf. In the 
early parts of this paper we presented a strategy to over- 
come these complications. 

The idealized runs we presented in section 7 have 
most in common with the simulations discussed by 
Grosfeld et al. (1997). Their setup differed slightly in 
that they imposed a w'ind field on the open ocean, and 
their initial conditions were an unstratified ocean with 
a temperature slightly colder than the one we have used. 
There are qualitative similarities between the results, 
namely, the pyclonic flow beneath the ice shelf and the 
general pattern of melting and freezing, but also distinct 
differences, particularly in the strength of the circula- 
tion. A cursory comparison with observation may sug- 
gest that the stronger circulation is closer to reality, but 
in nature there is considerable extra forcing on the sys- 
tem. For example, the growth and decay of sea ice drives 
a seasonal cycle of mixing and restratification in the 
water column to the north of the ice front, and the impact 
of this on flow beneath an ice shelf has been documented 
(Nicholls and Makinson 1998). Any process that excites 
motion in the water could promote more rapid vertical 
heat transfer, which would provide further forcing on 
the flow. Thus, the weak response to our idealized forc- 
ing may be correct. 

If the differences between our results and those of 
the earlier three-dimensional modeling studies lie in the 
different parameterizations of vertical heat transfer, we 
may conclude that a good representation of this process 
is essential for a realistic simulation of sub-ice shelf 
circulation. Our other principal conclusion is that the 
ice front presents less of a dynamic barrier to the buoy- 
ancy-driven circulation than has been suggested by the 
earlier studies. These differences highlight the need for 
thorough testing and intercomparison of a range of nu- 
merical models before more definitive statements can 
be made about the key processes controlling the cir- 
culation beneath ice shelves. We hope that the model 
we have presented here, both in its current form and 
with future developments, will prove a useful tool in 
advancing our currently limited understanding of these 
processes. 
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ABSTRACT 

Upper boundary conditions for numerical models of the ocean are conventionally formulated under the premise 
that the boundary is a material surface. In the presence of an ice cover, such an assumption can lead to 
nonconservative equations for temperature, salinity, and other tracers. The problem arises because conditions 
at the ice-ocean interface differ from those in the water beneath. Advection of water with interfacial properties 
into the interior of the ocean therefore constitutes a tracer flux, neglect of which induces a drift in concentration 
that is most rapid for those tracers having the lowest diffusivities. If tracers are to be correctly conserved, either 
the kinematic boundary condition must explicitly allow advection across the interface, or the flux boundary 
condition must parameterize the effects of both vertical advection and diffusion in the boundary layer. In practice, 
the latter alternative is often implemented, although this is rarely done for all tracers. 


1. Introduction 

Exchanges of heat and freshwater at the surface of 
the ocean represent the major forcing behind the crea- 
tion of new water masses. In the polar oceans, the source 
regions for most of the deep and bottom waters, the 
processes of exchange are strongly influenced by the 
presence of ice. If ocean models are to simulate the 
characteristics and distribution of abyssal waters cor- 
rectly, they must include an adequate representation of 
the coupling between ice and ocean. Because melting 
and freezing entail a transfer of water between the two 
media, this coupling formally manifests itself in the 
boundary conditions placed on both the vertical velocity 
and on the tracer concentrations. However, it is an al- 
most universal practice in numerical modeling of the 
ocean to regard the sea surface as a permeable, but 


Corresponding author address: Dr. Adrian Jenkins, British Ant- 
arctic Survey, Natural Environment Research Council, High Cross, 
Madingley Rd., Cambridge CB3 0ET, United Kingdom. 

E-mail: ajenkins@bas.ac.uk 


material, interface — properties may diffuse across it, but 
no water may be exchanged between the ocean and the 
overlying medium, be it atmosphere or ice. Huang 
(1993) criticized such practice, on the grounds that it 
suppresses the Goldsburgh-Stommel circulation and 
recommended setting the vertical velocity relative to the 
upper boundary equal to the precipitation, evaporation, 
melt, or freeze rate. Here we examine what the assumed 
nature of the ocean surface implies for the balance equa- 
tions for tracers under conditions of melting and freez- 
ing- . ■ L 

The ice-ocean boundary layer is a region that is dis- 
tinguished by relatively high gradients in tracer con- 
centrations combined with large changes in the effective 
diffusivity. It is often the case that the equations used 
to describe the transport of tracers in the bulk of the 
oceanic water column are not applicable in the boundary 
region. This problem can be circumvented by incor- 
porating a simple parameterization of the boundary lay- 
er processes into the boundary conditions. Such a for- 
mulation invariably entails the diagnosis of conditions 
at the ice-ocean interface that are distinct from those 
pertaining at the uppermost model grid point within the 
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sub-ice water column. In the presence of such differ- 
ences, the advection of meltwater into the bulk of the 
ocean plays an important role in the balance equations 
for tracers, but the advective flux is formally absent it 
the ice-ocean interface is treated as a material surface. 

We use a simple, one-dimensional model ot the ocean 
beneath sea ice to demonstrate the role played by melt- 
water advection in the conservation of tracers. To this 
model we apply two sets of boundary conditions that 
differ only in the assumption made about the nature of 
the ice-ocean interface: material or not. We also discuss 
the two-dimensional model of thermohaline circulation 
beneath an ice shelf presented by Hellmer and Olbers 
(1989) and Hellmer et al. (1998). This model has a 
material interface at the ice-ocean boundary, but chang- 
ing the kinematic boundary condition would entail a 
considerable increase in complexity. Instead, we include 
the effect of meltwater advection through a modified 
flux boundary condition and show' that this simple fix 
improves the model performance. Finally, we discuss 
how the boundary formulations discussed here compare 
with those used elsewhere. 

2. Fundamental equations 

Our discussion will focus on the evolution of scalar 
properties in the upper ocean, and we will use the fol- 
lowing generic equation to model this evolution: 

— + V • (uX) + — (wX) = -v • (u'X') - f (w'X'), 
dt dZ °Z 

( 1 ) 

where X could be temperature T f salinity S, or any other 
conservative tracer; boldface variables represent hori- 
zontal vectors; primed variables represent turbulent fluc- 
tuations; and angle brackets indicate ensemble averag- 
ing. In the presence of an ice cover, this equation is 
forced by the upper boundary conditions placed on the 
terms describing the vertical advection and diffusion of 
tracers. 

The advective flux at the ice-ocean interface is simply 


W* = d jf + U„ • Vz*. (4) 

Under this assumption the melt/freeze rate does not im- 
pact the continuity equation, a simplification that ex- 
plains why (4) is generally favored over (3). The price 
to be paid for this simplification is a slight inconsistency 
between the continuity and tracer equations because the 
impact of melting and freezing must be retained in the 
latter. 

We express the diffusive flux at the ice-ocean inter- 
face as the product of a turbulent transfer coefficient 
and the difference in the tracer concentration across the 
surface boundary layer: 

<wX)L - y x (X - xj, (5) 

where the transfer coefficient y x has the dimensions of 
velocity (i.e., diffusivity over length). Since the surface 
boundary layer itself is part of the ocean domain, the 
precise meaning of the transfer coefficient depends on 
the form of discretization applied to (1). In a layer mod- 
el, X in (5) will be a depth-averaged value across the 
upper layer, part of which will be the boundary layer, 
and the value of y x should reflect this. In a level model, 
X will be the value at the first level below the boundary, 
and the value of y x should reflect whether this point is 
within or outside of the boundary layer. At vertical res- 
olutions commonly used in ocean models, these dis- 
tinctions are unimportant because most of the change 
in tracer concentration occurs across a very thin viscous 
sublayer (McPhee 1990). The result is that conditions 
at the upper level, or within the upper layer, of a model 
are practically equivalent to the far-field conditions out- 
side the boundary layer. 

The application of Eqs. (2) and (5) requires knowl- 
edge of the tracer concentration at the ice-ocean inter- 
face. We derive this from a consideration of the budget 
within an infinitesimal control volume that straddles the 
interface (Mellor et al. 1986). We demand an exact bal- 
ance between the vertical diffusive fluxes and the sourc- 
es or sinks of the tracer represented by the consumption 
or production of meltwater within the control volume: 


(wX)\ b = w b X b , (2) 

where the subscript b indicates conditions at the ice- 
water contact. For a nonmaterial surface, the vertical 
velocity at the boundary is derived from a combination 
of the motion of the boundary and the advection of fluid 
across the boundary: 

W b = ^ + n b ■ \z h - m, (3) 

ot 


- P7x(X h “ X) = pm{X b - X,). (6) 

Here the subscript i indicates ice properties, p are den- 
sities, and & is a diffusivity. For most tracers, including 
salinity, stable isotopes and incorporated gases, the dif- 
fusivity in ice is negligible, so the first term on the left- 
hand side of (6) vanishes. This is not so in the case of 
temperature, for which Eq. (6) becomes 



where m indicates the melt rate (negative for freezing) 
expressed as a thickness of seawater per unit time. For 
the more commonly encountered material interface, the 
particles lying on the boundary at one time must lie 
there at all subsequent times, so their motion must be 
equal to that of the boundary itself: 



pcy r (T b ~ T) - pmL it 


(7) 


where c are specific heat capacities and L a latent heat 
of fusion. 

Complete specification of the boundary problem re- 
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Fig. 1. Schematic of a one-layer ocean beneath sea ice. The ocean 
layer is characterized by a single temperature and salinity ( T , S) and 
is forced by an atmospheric heat flux (Q h ) that passes through leads. 
The temperature and salinity at the ice-ocean interface ( T b , S h ) differ 
from those in the interior, and the differences drive melting and freez- 
ing. The sea ice is regarded as a continuum, so the ice-ocean bound- 
ary conditions and the atmospheric heat flux are formally applied at 
the level of the mean ice draft ( z = -h). The seabed lies at z ” -H. 


quires an expression for the melt/freeze rate. This is 
determined by the rate at which heat and salt are sup- 
plied to, or extracted from, the control volume, and can 
be calculated from the combination of (7), with an ap- 
propriate specification or parameterization of heat con- 
duction through the ice, and the salinity version of (6). 
A third constraint is that the seawater within the control 
volume must be at the freezing point, which we define 
using a linear equation: 

T b = A t S b + A 2 + A 3 P„ (8) 

where P b is the pressure at the interface and A,_ 3 are 
constants: -0.0573°C/psu, 0.0832°C, and -7.53 X 10 8 
°C/Pa, respectively. This relationship is based on that 
of Millero (1978), with the salinity dependence line- 
arized. Equations (6)-(8) can be solved simultaneously 
for m, S b , and T h , a procedure that is simplified by the 
use of a linear equation in (8). With the melt rate known 
the boundary concentration of any other tracer can be 
calculated from (6) if the concentration in the ice is 
known. 

Boundary conditions of the form given in (5)-(8) 
were introduced by Josberger (1983) and McPhee 
(1983). Mellor et al. (1986) recognized the importance 
of molecular diffusion within the viscous sublayer in 
determining the melt/freeze rate, and McPhee et al. 
(1987) and Steele et al. (1989) investigated expressions 
for the turbulent transfer coefficients that included ex- 
plicit parameterizations of this process. The magnitudes 
of the transfer coefficients used in this study are derived 
from the results of this pioneering work. 

3. A simple model of the ocean beneath sea ice 

We begin by considering the evolution of a one-layer 
ocean model (Fig. 1) subject to a symmetric cycle of 


heating and cooling. We assume horizontal homogeneity 
and treat the lower boundary as a fixed, insulating sur- 
face, thus isolating the upper boundary condition as the 
only forcing on the model. Under these assumptions 
integration of Eq. ( 1 ) over the depth of the water column 
leads to 

’* — dz + w b x„ = -<»'*%. (9) 

H &t 

The integral expression can be evaluated as 

f —dz = j[(H - /i)X] + X b j (10) 
dt dt dt 

where the overbar indicates a depth-averaged quantity, 
and we have made use of the Leibnitz theorem for the 
differntiation of an integral. 

If there were complete ice coverage, the turbulent 
fluxes at the upper boundary would take the form given 
in (5) but, if direct ocean-atmosphere exchange can take 
place through leads, the flux boundary condition at the 
ice-ocean interface (z = —h) becomes 

piw'XX = Apy x (X - X,) + (1 - A)Q X , 0D 

where A is the areal fraction of the surface covered by 
ice and Q x is the tracer flux passing through the leads. 
In the temperature equation, Q x is replaced by QJc , 
where Q h is the heat flux. The two choices for the re- 
maining advective boundary condition are discussed 
separately below. 

We force the model with a sinusoidal variation in the 
atmospheric heat flux having an amplitude of 500 W 
mr 2 over open water and a period of one year. For all 
other tracers we assume that the flux passing through 
the leads is zero. The water column thickness (H - h) 
is initially 50 m and the ocean has an initial salinity of 
34.5 psu and an initial temperature at the appropriate 
salinity-dependent freezing point. The initial ice thick- 
ness is arbitrary, as we are only concerned with how 
the melting and freezing cycle influences the bulk ocean 
properties. The areal coverage of ice is assumed constant 
at 90% and the model is run for 10 years. We use a 
turbulent heat transfer coefficient, y r , of 5 X 10 s m 
s" 1 and a turbulent salt transfer coefficient, y v , of 
0.04 y r . These values are appropriate for a moderate 
friction velocity of 5 X 10 3 m s' 1 (McPhee et al. 1987; 
Holland and Jenkins 1999). 

For the purposes of diagnosing the melt/freeze rate 
and the temperature and salinity at the ice-ocean bound- 
ary, we assume that no salt is incorporated into the ice 
on freezing. We also assume the ice to be sufficiently 
thick that the heat conducted through it is negligible, 
which allows us to set the first term on the left-hand 
side of (7) to zero. The pressure-dependent term in (8) 
can also be neglected. Note that these assumptions, 
along with the simple formulations we have used for 
the heat and salt transfer coefficients, do not affect the 
form of the equations presented below. They influence 
only the magnitude of the forcing. 


v 
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a. With meltwater advection 

We first regard the ice-ocean boundary as a non- 
malerial interface and apply a kinematic boundary con- 
dition of the form given in (3): 

w* = “ “ Am - (12) 

Applying this and the flux boundary condition given in 
( 1 1 ) to (9), the conservation equation for a generic tracer 
becomes 

n — Qx 

— [(H - /i)X] - AmX h = Ay x (X h - X) - (1 - A)—. 

(13) 

The vertical velocity must vanish at the seabed, so our 
assumption of zero gradient in the horizontal velocity 
field implies that the w b must also be zero, and (12) can 
be rewritten: 


—(H - h) = Am (14) 

dt 

because H is fixed for all time. Using (14) to expand 
the derivative on the left-hand side of (13), the model 
equations for salinity and temperature, respectively, be- 
come 


dS _ m ) 

Tt ~ (H - h) 


S) 


(15a) 


dT _ A(y T + m) 
dt ~ (H - h) ( h 


f - 11 ~ Mh 

’ pc(H - h y 


(15b) 


Results obtained with this version of the model are 
shown in Fig. 2. The annual cycle in ice thickness shows 
a peak-to-peak range of about 1 .5 m, giving rise to ocean 
salinities that cycle over a range of about 1 psu. Ocean 
temperatures deviate by up to 0.35°C from the freezing 
point. The maxima and minima in ice thickness occur 
about 2.5 weeks after the change in sign of the atmo- 
spheric heat flux in spring and autumn (Fig. 2a). The 
lag is approximately equal to the reaction time, t, of 
the oceanic layer: 


(H - h) 

T = , 

TVefr 


(16) 


where y /r)1 is the effective heat transfer coefficient, de- 
fined by Holland and Jenkins (1999) as 


l . A,(S t - S,)cy , 

y ’« ‘ H Wr, 

Continual relaxation towards the freezing point means 
that the water temperature is a function of the rate at 
which heat is lost to, or gained from, the atmosphere, 
subject to the lag mentioned above, so the temperature 
extremes occur shortly after midsummer and midwinter 
(Fig. 2b). The salinity is a function of the cumulative 





Fig. 2. Results of a ten-year integration of (15): (a) change in ice 
thickness, (b) temperature, and salinity of the ocean. In (a) the dotted 
lines indicate the time at which the atmospheric heat flux changes 
sign in spring (S) and autumn (A). In (b) the stars labeled S (spring), 
S (summer), A (autumn), W (winter) indicate the times of maximum 
and zero atmospheric forcing, while the straight line labelled T f in- 
dicates the freezing point at atmospheric pressure. 


heat loss or gain throughout a growth or melt season, 
so the extremes occur when freezing switches to melt- 
ing, and vice versa. Together these signatures give an 
annual cycle that traces an ellipse (counterclockwise) in 
temperature/salinity space. This cycle does not exactly 
repeat itself. There is a drift of about 0,0025 psu yr 1 
toward lower salinities, consistent with a reduction in 
ice thickness of approximately 4 mm yr -1 . This drift 
arises because the ice cover is thicker on average during 
the spring to autumn period of heating than during the 
autumn to spring period of cooling (Fig. 2a). Tn this 
model formulation the water lost by the ice during melt- 
ing is gained by the ocean, and vice versa during freez- 
ing. The water column is therefore thinner and has a 
smaller heat capacity, on average, during periods of 
heating than during periods of cooling, and this allows 
slightly more melting to occur than freezing. 

Despite the simplicity of the model we find some 
support for such a cycle of ocean properties in obser- 
vations made during the Arctic Ice Dynamics Joint Ex- 
periment. Figure 3 shows mixed layer temperatures and 
salinities recorded daily between late May 1975 and 
mid- April 1976 at the Snowbird station (Maykut and 
McPhee 1995). Many factors not considered in our one- 
dimensional model, such as drift of the station, hori- 
zontal advection of water masses, mixing across the 
pycnocline, variatlQlIs in ice concentration, and a var- 
iable input of turbulent kinetic energy to the ocean with 
the passage of weather systems, will all have influenced 
the observations. Given these complicating factors, the 
basic shape and phasing of the observed annual cycle 
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Fig. 3. A 331 -day record of temperature and salinity at the Snow- 
bird ice camp. Each dot represents an average value for the upper 
25 m of the water column. Data from the first day of each month are 
circled and labeled. The straight line labeled 7} indicates the freezing 
point at atmospheric pressure, and the two dashed lines indicate the 
likely uncertainty in the observed deviation from the freezing point 
(Maykut and McPhee 1995). 


in mixed layer properties compare favorably with the 
model results. 

The main qualitative difference between the obser- 
vations and model results is the lack of any sustained 
period of supercooling, at levels above the detection 
threshold of the instrumentation used, in the observa- 
tional record. A lack of supercooling suggests that the 
ocean response time is shorter during periods of ice 
growth than during periods of melting. The reason may 
be frazil ice production in the water column or funda- 
mentally different behavior within the boundary layer 
when salt, rather than freshwater, is produced at the 
phase change interface. Whatever the physical mecha- 
nism, we can simulate the effect in the model by in- 
creasing the magnitude of the heat and salt transfer co- 
efficients during periods of supercooling. The results 
are shown in Fig. 4. With the wintertime response of 
the ocean 20 times faster than before, supercooling is 
restricted to no more than about 0.02°C, and the peak 
in supercooling and the transition from melting to freez- 
ing are practically coincident with the maximum and 
the change in sign, respectively, of the atmospheric forc- 
ing. Because there is little phase lag between the cycles 
of atmospheric heat flux and ice growth during half of 
the year, the drifts in ice thickness and ocean salinity 
are reduced by a factor of 2. 

b. Without meltwater advection 

We next treat the ice-ocean boundary as a material 
surface and apply a kinematic boundary condition anal- 
ogous to (4): 




Fig. 4. Same as Fig. 2 but with heat and salt transfer coefficients 
increased by a factor of 20 when the ocean temperature falls below 
the freezing point. 




dh 

dt' 


(18) 


Applying this and the flux boundary condition given in 
(11) to (9) gives 

" h)X] = Ay x {X h - X) - (1 - A)— (19) 

dt P 

for a genetic tracer. Equating w b to zero yields 


-(H - h) = 0, 
dt 


( 20 ) 


and the model equations for salinity and temperature 
become 


dS 

dt 


Ay s 

D 


(S h ~ S) 


(21a) 


K - — T\ — (! - VQ* 

dt D { b } pcD 


(21b) 


where D = H - h is the, now temporally constant, 
water column thickness. 

The main differences between (15) and (21) arise 
because the water advected into or out of the ocean in 
(15), which is absent in (19), has properties equal to 
those at the boundary. Advection therefore constitutes 
a tracer flux that is additional to the diffusive fluxes 
given by (5). Using (21) we obtain a similar evolution 
of the ice cover to that obtained with (15), but we find 
a rapid drift of about 0.1 psu yr _! toward higher salinity 
(Fig. 5). The increase in salinity is what we would an- 
ticipate for a net ice growth of 15 cm yr“ \ and is clearly 
in error. It is a result of the missing advective fluxes, 
which are expressed as the products of the melt rate and 
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Fig. 5. Results of a ten-year integration of (21): (a) change in ice 
thickness, (b) temperature, and salinity of the ocean. In (a), the dotted 
lines indicate the time at which the atmospheric heat flux changes 
sign in spring (S) and autumn (A). In (b), the stars labeled S (spring), 
S (summer), A (autumn), W (winter) indicate the times of maximum 
and zero atmospheric forcing, while the straight line labeled T } in- 
dicates the freezing point at atmospheric pressure. 

the difference in tracer concentration between the mixed 
layer and the ice-ocean interface [Eq. (15)]. The con- 
centration differences are generated by the rejection or 
uptake of tracers during freezing or melting, so they 
change sign synchronously with the melt rate (Fig. 6). 
Hence, the advective fluxes are always of the same sign, 
and, although they are generally small, their neglect 
leads to an error that accumulates throughout the annual 
cycle. With the reaction time of the ocean shortened by 



Fig. 6. Illustration of the net dilution induced by advection of 
meltwater into an ocean layer having a nonmaterial upper boundary: 
(a) melt rate, (b) salinity difference across the ice-ocean boundary 
layer, and (c) salinity change that results from the inclusion of melt- 
water advection on the right-hand side of (15a), defined by dS/dt - 
Am(H - h)~'(S k - S ). Dotted lines indicate the times at which the 
melt rate changes sign. 



Time (yr) 



Salinity (psu) 

Fig. 7. Same as Fig. 5 but with heat and salt transfer coefficients 
increased by a factor of 20 when the ocean temperature falls below 
the freezing point. 

a factor of 20 during freezing, the difference in salinity 
across the boundary layer, and hence the salt flux error, 
is reduced in winter. The annual drift in salinity is there- 
fore approximately halved (Fig, 7). 

Although the net salinity drift is a function of our 
model setup, in particular the choice of water column 
thickness, the magnitude of the flux errors are deter- 
mined only by the nature of the ice-ocean boundary 
conditions. In the appendix we develop general expres- 
sions for the size of the errors and demonstrate that they 
are approximately proportional to the square of the de- 
viation in water temperature from the freezing point. 
The errors are most significant for slowly diffusing trac- 
ers with small y x , the ratio of the excluded advective 
to included diffusive fluxes being ml y x . For this reason 
the salinity error is much more apparent in Figs. 5 and 
7 than the temperature error. The solution could not drift 
from the freezing point line in any case, but the heat 
flux error alone would act to reduce the thickness of the 
ice cover and freshen the mixed layer. 

4. A model of thermohaline circulation beneath an 

ice shelf 

We next consider the model of thermohaline circu- 
lation beneath an ice shelf described by Hellmer and 
Oibers (1989). The model domain consists of a vertical 
plane, oriented perpendicular to the ice front and ex- 
tending into the cavity beneath the ice shelf. The bound- 
ary located at the ice front is open, while the interior 
boundaries are fixed material surfaces, the upper one 
corresponding to the ice shelf base. A single equation 
for the vertical streamfunction replaces the momentum 
and continuity equations, while the generic equation for 
tracers, equivalent to our Eq. (1), is written 
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^ + |-(uX) + |-(vvX) 
dt dy dz 


— I I + ~l Kv 


dy\ dy I dz\ dz 


dX 


+ C, 

(22) 


where v and vv are horizontal and vertical components 
of velocity, K H and K v are horizontal and vertical eddy 
diffusivities, and C indicates a convective term used to 
eliminate static instability. Boundary conditions on the 
streamfunction are specified such that there is no normal 
flow at the solid boundaries and only normal flow at 
the open boundary, where the model is forced with an 
observed vertical profile of temperature and salinity. 
Where outflow occurs, the observations are overwritten 
by values advected from the cavity interior. Heat and 
salt balances at the ice shelf base are specified as in 
Eqs. (6) and (7), with the first term on the left-hand side 
of (7) parameterized in a manner appropriate for heat 
conduction into a melting ice shelf (N0st and Foldvik 
1994; Holland and Jenkins 1999): 

= pmcXT , - T h l (23) 

where T s is the temperature at the surface of the ice 
shelf. 

The model equations are solved on a rectangular grid 
with regular horizontal and vertical spacing. Staggering 
of the grid means that the vertical velocity is defined 
on points that lie half a grid cell above and below the 
points carrying the tracers. The ice-ocean boundary 
passes through the velocity points and, since the tracer 
fluxes (K v dX!dz) are naturally defined on these points, 
the diffusive flux boundary condition can simply be 
written 

= y x (X b - X U2 ), (24) 

where the subscript 1/2 indicates conditions at the up- 
permost grid point within the water column, half a grid 
box below' the interface. The lack of tracer points on 
the boundary is not a problem because the value X b is 
derived from (6) and need not be formally assigned to 
any grid point. 

Application of (24) is equivalent to (21) in that, con- 
sistent with the assumption of no normal flow at the 
boundary, it neglects the advection of meltwater across 
the ice-ocean interface. However, if the tracer budgets 
are to be closed correctly, the first grid point in from 
the boundary must feel the additional sinks of heat and 
salt associated with the advection of meltwater across 
the boundary layer. Because Eq. (22) is cast in flux form, 
it would be simple to impose an advective flux at the 
ice-ocean interface, but the advection field w'ould then 
be inconsistent with the velocity field, which is derived 
on the basis of no normal flow at the boundaries. The 
problem of nonconservation of tracers would be exac- 
erbated. Instead we apply a modified flux boundary con- 





Fig. 8. Potential temperature vs salinity diagram comparing ob- 
servations made at the calving front of Pine Island Glacier (dots) 
with the output from a two-dimensional model of thermohaline cir- 
culation in the sub-ice cavity (Hellmer et al. 1998). The vertical 
profiles of potential temperature and salinity used here constitute the 
open boundary of the model domain. Model results obtained using 
flux boundary conditions defined as in (24) and (25) are indicated 
by open circles and crosses, respectively. Both versions of the model 
were tuned to give net melt rates for the glacier that were consistent 
with observation. 


dition that is analogous in form to the terms appearing 
in our earlier layer- integrated equations 15: 

= (Yx + rn){X b - X J/2 ). (25) 

Hellmer et al. (1998) applied the two-dimensional 
model to the water column beneath the 90-km floating 
extension of Pine Island Glacier, a fast-flowing outlet 
glacier of the West Antarctic Ice Sheet. Employing flux 
boundary conditions of the form shown in (25), they 
were able to obtain good agreement between modeled 
outflow characteristics and observations made at the 
calving front of the glacier. Figure 8 shows a comparison 
between the results of this model run and the obser- 
vations. The results of a second model run, which em- 
ployed boundary conditions of the form shown in (24), 
are also shown in Fig. 8. These latter results are clearly 
inferior. 

The straight line trajectory in potential temperature/ 
salinity space, apparent in Fig. 8, is characteristic of ice 
melting in seawater. Gade (1979) calculated the gradient 
of this line to be 

dT = c(T w - T h ) + Lj + c t (T h - TJ 
dS c(S„ - S,) 

where the subscript vv indicates the far-field seawater 
properties and for glacier ice the salinity is zero. The 


dX 

K v- 

dz 
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slope is primarily a function of the water and ice prop- 
erties and depends only very weakly on the melt rate, 
through its influence on the boundary temperature. Tak- 
ing the ratio of heat flux to salt flux, as expressed by 

(25) , substituting from (6), (7), and (23), and setting the 
salinity of the ice to zero, an expression analogous to 

(26) is readily obtained. However, writing the fluxes as 
in (24), a similar procedure leads to 

d_T = h + C ‘ iT * - 12 (27) 

as c(S b - s,) ’ 

which represents erroneous forcing on the model. Fur- 
thermore, the appearance of the boundary salinity in the 
denominator makes the slope sensitive to such factors 
as the melt rate and the magnitude of the turbulent trans- 
fer coefficient for salt. 


5. Significance for other models 

Boundary formulations that include all the effects of 
meltwater advection are, to our knowledge, uncommon. 
However, for conditions of low melt/freeze rate and for 
integrations over short timescales, the errors caused by 
nonconservation of tracers are probably small. This is 
the case with most models of the interaction between 
ice shelves and the ocean, the early versions of which 
(Hellmer and Olbers 1 989; Scheduikat and Olbers 1 990; 
Jenkins 1991) all used boundary formulations analogous 
to (21) or (24). Later versions of these models (Nicholls 
and Jenkins 1993; Jenkins and Bombosch 1995; Hellmer 
et al. 1998) have been corrected, although it is only in 
the case of the Pine Island Glacier simulations cited in 
the previous section where melt rates are ~~10 m yr \ 
an order of magnitude higher than in the other examples, 
that the correction has had a significant impact on the 
results. 

The literature on coupled sea ice-ocean models is too 
vast to enable us to present an authoritative survey of 
all the various boundary formulations in use. In any 
case, the published descriptions of the models often lack 
the level of detail required to reveal the use of conser- 
vative or nonconservative boundary conditions. Our dis- 
cussion has centered on what Holland and Jenkins 
(1999) refer to as three-equation formulations. These 
employ all of Eqs. (6)-(8) to diagnose the temperature, 
salinity, and melUfreeze rate at the ice-ocean interface. 
This type of formulation is not commonly encountered 
in sea ice models, but where it is, boundary conditions 
corresponding to (21) and (24) appear to be in use. Two- 
equation formulations are encountered more often. In 
these, the salinity at the boundary is assumed equal to 
that at the upper model level or layer, and only heat 
transfer through the thermal boundary layer is consid- 
ered. McPhee et al. (1999) demonstrate that with a suit- 
able choice of effective heat transfer coefficient, such a 
formulation can yield realistic melt rates. Since the 
boundary salinity is not explicitly calculated, the con- 


ventional choice for the salinity equation is a flux 
boundary condition of the form 

= m(S, - 5), (28) 

where S represents the salinity at the uppermost grid 
point in the water column. This is simply the expression 
in (25), rewritten using (6), so it is a conservative bound- 
ary condition that implicitly includes both advective and 
diffusive fluxes. However, the analogous boundary con- 
dition for temperature, obtained from a combination of 
(25) and (7), 



Kv 


dT 

dz 


b 



(29) 


is rarely seen. The conventional version of (29) includes 
only the latent heat term and the heat conduction through 
the ice, and therefore conforms to the nonconservative 
form of (24). Incorporation of meltwater advection into 
the salt balance eliminates the dominant source of error, 
but the advection term can rise to several percent of the 
latent heat term in (29). Although uncertainties in the 
latent heat term itself may still dominate the overall error 
in the heat flux, errors resulting from such uncertainty 
are of opposite sign for melting and freezing, whereas 
the error arising from neglect of the advective term is 
independent of the sign of the phase change. 

More generally, ice-ocean models often include ad- 
ditional processes such as percolation of surface melt- 
water through to the ice-ocean interface and precipi- 
tation and evaporation at the surface of leads. If per- 
colating meltwater is assumed to have the salinity of 
the ice and be at the freezing point, Eqs. (6) and (7) 
become 

py s (S - S b ) = p(m + m p )(S b - S t ) (30) 
+ pcy T (T- T b ) = pmL, 

+ pcm p (T b - 7},), (31) 

where m p is the percolation velocity, which is either 
specified or derived from an ice model, and T fi is the 
freezing point appropriate for a salinity of S t . The ver- 
tical velocity at the ice-ocean interface, taken as the 
level of the mean ice draft, is expressed as 


p,c t k ( 


dj\ 

dz 


fob , 

— + IK 

dt b 


Vz h - A(m + m ) 


- (1 - A)(p - e), (32) 

where p and e represent the rates of precipitation and 
evaporation respectively. We could now either apply 
separate advective and diffusive boundary conditions 
based on these, as in section 3, or, if we wished to 
maintain the simplicity of a material ice-ocean inter- 
face, we could apply a combined flux boundary con- 
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dition of the form given in (25). For salinity this would 
read 

= A(y s + m + - S) 

- (1 - A)(p - e)S. (33) 

Again we emphasize, by substitution from (30), that this 
yields the conservative boundary condition: 

= A(m + m,)^. - 5) - ( 1 - A)(p - e)S. (34) 

b 

In contrast, neglect of the meltwater advection terms in 
(33) yields a boundary condition analogous to (24), 

= Ay s (S b - S) - (1 - A){p - e)S, (35) 

b 

and leads, on substitution from (30), to the nonconser- 
vative form: 

= A(m 4- m p ){S t - S h ) - (1 - A){p - e)S. (36) 

b 

In a two-equation formulation (34) would be the natural 
choice, but we note that the advection of percolating 
meltwater into the ocean has an impact on the ice-ocean 
heat flux, equivalent to the addition of the term m p {T fi 
- T) to the right-hand side of (29), that is likely to be 
overlooked. This term would generally be small but, 
since the interfacial melt rate and percolation velocity 
are unrelated, it is possible that under certain conditions 
it could be the dominant term in the ice-ocean heat flux. 

Larger-scale ice-ocean models often use a simpler 
one-equation formulation in which only the freezing 
point at the upper level or layer of the model is diag- 
nosed. At each time step the computed temperature is 
reset to the freezing point and an appropriate amount 
of ice melted or frozen. Once again, because a separate 
boundary salinity is not calculated, the salinity balance 
is almost certainly computed correctly, but we suspect 
that the much smaller errors in the heat balance are 
routinely made. Of course, none of the errors discussed 
in this paper will be apparent in models that employ 
either surface restoring fluxes or diagnosed flux correc- 
tions, which are specifically designed to eliminate drift 
such as that in Figs. 5 and 7. 

6. Closing remarks 

Mellor et al. (1986) concluded that for small |m/w*|, 
where n* is the surface friction velocity, it is possible 
to neglect the influence of vertical advection in the for- 
mulation of boundary conditions at an ice-ocean inter- 
face. However, this conclusion was based on an as- 
sessment of the likely impact of vertical advection on 
the structure of the boundary layer. We do not dispute 
this finding and have implicitly made use of it in our 






parameterization of diffusion through the boundary lay- 
er. Nevertheless, we conclude that meltwater advection 
does play a role in the balance equations for tracers 
within a sub-ice water column and that this role may 
be neglected only if \mfy x \ is small. Since y 5 /w* ** 4 
X 10 4 and the turbulent transfer coefficients for most 
other dissolved species are of a similar order of mag- 
nitude to that of salt, this latter criterion is a much more 
stringent requirement. We conclude that in general it is 
not possible to neglect the impact of meltwater advec- 
tion. 

We have presented a formulation of the flux boundary 
condition at an ice-ocean interface that includes both 
advective and diffusive fluxes and have demonstrated 
that inclusion of the advective term is necessary for the 
conservation of tracers. In practice, the advective fluxes 
are likely to be excluded only if the application of the 
boundary conditions involves the explicit diagnosis of 
the tracer concentration at the ice-ocean interface. There 
are relatively few models that include this level of so- 
phistication and, although the heat flux error is an ex- 
ception in that its occurrence is likely in almost all for- 
mulations, it is also the least significant of the errors. 
For other tracers, the commonly used boundary con- 
ditions already incorporate both advective and diffusive 
fluxes, but we believe the foregoing discussion to be of 
more than academic interest. As computing resources 
become ever more readily available, more detailed rep- 
resentations of key physical processes are likely to be- 
come part of even global-scale models. While our cur- 
rent knowledge of ice-ocean boundary physics is not 
sufficient to say with confidence that melting and freez- 
ing can be calculated more accurately if the boundary 
salinity is explicitly diagnosed, the inclusion of some 
additional tracers, in particular stable isotopes, does re- 
quire the computation of the their boundary concentra- 
tions. Also, if the use of a nonmaterial interface at the 
ocean surface, as recommended by Huang (1993), be- 
comes more commonplace, the distinction between the 
advective and diffusive parts of the tracer fluxes will 
become important. If a separate advective boundary 
condition is applied to the scalar equations, the purely 
diffusive flux boundary condition (24) must be em- 
ployed to avoid double-counting of the meltwater ad- 
vection term. Finally, as efforts are already being made 
to eliminate the need for flux corrections in coupled 
ocean-atmosphere general circulation models, it is im- 
portant to isolate all possible causes of spurious drift in 
the properties of the ocean. 
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APPENDIX 

How Large Are the Errors? 

The error in the tracer flux, per unit area of ice cover, 
associated with the use of boundary conditions of the 
form given in (21) and (24) is 

Q xm = -pm(X b - X), (Al) 

which, using (6) and ignoring the diffusive flux into the 
ice, can be rewritten: 

Qx~ = P — <*» - X'l (A2) 

For tracers that are rejected from the ice on freezing, 
the drift will be toward higher concentrations. For 
equivalent melt/freeze rates the drift will tend to be 
slightly faster during freezing, when the boundary con- 
centration is raised, than during melting, when the 
boundary is depleted in the tracer. By substituting for 
the melt rate from (7) we can express the error as 


y T = 0,«*. ( A6 ) 

and setting the ratio y s f y T to a constant value. We take 
values of 0.0 1 for the factor C h and 0.04 for the ratio 
y s /y T (McPhee et al. 1999; Holland and Jenkins 1999). 

We consider first a situation, such as that addressed 
in the model described in section 3, in which the heat 
conduction term is negligible and to a good approxi- 
mation we can assume melting and freezing to be driven 
entirely by temperature variations in the water column. 
Under these conditions, we can ignore the first terms 
on the left-hand side of both (6) and (7) and hence obtain 
a simple expression for the ratio of the salinity to tem- 
perature difference across the boundary layer: 

(S b ~ S) = (. S b - S t )cy T 
(T b - T) L t y s 

We refer to the elevation of the far-field temperature 
above the freezing point as the thermal driving, 

T* = T - A,5 - A 2 - A 3 P b , (A8) 

and, using (8), we relate this to the temperature differ- 
ence across the boundary layer: 


T„ - r= a ,(s, - S) - r*. 

(A9) 

Using (A7) we can rewrite this as 


1 

II 

(A 10) 

where 


. . A,(5 fc - S,)cy T 

A = “ L, y, ■ 

(All) 


This rearrangement is the origin of the effective heat 
transfer coefficient, defined in (17). With (A6) and 
(A 10) substituted into (A3) and (A5) we obtain 


Qxctt 


(X b - X,) 

PJxV 


pcy r (T - T b ) + p i c i 



(A3) 


Qxc IT 


C h (X b -X t ) ( pc „ 

w* — - / * 

piyxhr) \ L A 


2 


and 


(A 12) 


The error in the heat flux is expressed 

Qhcn = ~pcm(T b - T), (A4) 

and, substituting directly for the melt rate, we obtain 


The external factors that determine the melt/freeze 
rate are the level of turbulence in the mixed layer, how 
far the temperature of the mixed layer deviates from the 
freezing point, and the temperature gradient in the ice. 
We focus on the salt and heat flux errors and seek to 
express (A3) for salinity and (A5) in terms of the above 
three factors. The level of turbulence determines the 
magnitude of the transfer coefficients for heat and salt, 
an effect that we parameterize by making y T a linear 
function of the friction velocity, 



where the temperature gradient terms within the paren- 
theses of (A3) and (A5) have been set to zero, consistent 
with our assumption of negligible heat conduction. 

Holland and Jenkins (1999) demonstrate that for ther- 
mal driving of <0.5°C, the small variations in boundary 
salinity have little impact on the effective heat transfer 
[Eq. (17)]. Taking S b in (All) to be constant at 34.5 
psu, and taking the ice salinity to be zero, A has a 
constant value of 1.59. The only remaining variables in 
(A 12) and (A 13) are then the thermal driving, the fric- 
tion velocity, and the difference in tracer concentration 
between the ice and the boundary. To evaluate the salt 
flux error we make the same approximations for the 
boundary and ice salinities. We then replace all con- 
stants in (A 12) and (A 13) with their numerical values 


pcy r (T - T b ) + PiC t kr~ 
oZ 


■ (A5) 
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x 10~ 3 



Fig. Al. Salt and heat flux errors induced by the neglect of meltwater advection across an ice- 
ocean boundary for friction velocities of 0.1, 0.5, 1.0, and 2.0 cm s~': (a) salt flux errors defined 
by (A 12) (solid line) and by the approximation (dashed line) given in (Al4a), (b) heat flux errors 
defined by (A 13) (solid line) and by the approximation (dashed line) given in (A 1 4b), (c) salt 
flux errors defined by (A 17) (solid line) and by the approximation (dashed line) given in (A 19a), 
and (d) heat flux errors defined by (A 18) (solid line) and by the approximation (dashed line) given 
in (A19b). 


and express the salt and heat flux errors, respectively, 
as 

Q Sc „ * 0 .5u*Tl [psu kg s~' nr 2 ] (A14a) 
Q hc „~ 200u*Tl [Wm- 2 ], (A14b) 

where h* is in meters per second and 7* is in degrees 
Celsius. 

We next consider the alternative situation where the 
thermal driving is zero and freezing proceeds as a result 
of the heat conducted into the cold ice above. Under 
these conditions we can combine (6) for salinity with 
(A9) to give an expression for the melt rate, 


m 


ys 

A AS* - S') 


(T ~ T b \ 


(A 15) 


and substituting this into (7) we can express the tem- 
perature difference across the oceanic boundary layer 
in terms of the temperature gradient at the base of the 
ice: 


T~T b 


(1 - A) dT t 
A pcy T dz b 


(A 16) 


With (A6) and (A 16) substituted into (A3) and (A5) we 
obtain 


(X„-X.) 1 Uc&dT, y 

pC,i(y x /y r ) \ A dz h J 


and (A 17) 


Qhcxr 


(l - A) 1 ( P'C'k'dT ; 


pC h L { w* l A dz 


(A 18) 


Making the same approximations as before, the salt and 
heat flux errors can be written, respectively 


[psu kg s" 1 m 2 ] (A19a) 

[Witt 2 ], (A 19b) 


1.3 x 10 ’i 

<dT, 

W* * 

K' iz 

-3 X 10 2 

(<>T, 

w* 

U 


where, once again, u * is in meters per second and the 
temperature gradient is in degrees Celsius. 

The salt and heat flux errors given by (A 12), (A 13), 
(A 14), (A 17), (A 18), and (A19) are plotted in Fig. Al 
for a range of values of friction velocity, thermal driv- 
ing, and basal temperature gradient. The quadratic de- 
pendence on the latter two factors means that the errors 
are always of the same sign. In the case of (A 14) the 
effects of the errors are counteractive, in that the heat 
flux error should cause melting and dilution of the mixed 





296 


JOURNAL OF PHYSICAL OCEANOGRAPHY 


Volume 31 


layer, while the salt flux error acts to increase the mixed 
layer salinity. In the case of (A 19), both errors serve to 
accentuate the effects of freezing. We note that however 
the errors combine and whichever is dominant, the effect 
is most evident in the computed salinity because the 
properties of the seawater are constrained to follow 
closely the liquidus, which has an almost flat trajectory 
in temperature/salinity space. 

REFERENCES 

Gade, H. G., 1979: Melting of ice in sea water: A primitive model 
with application to the Antarctic ice shelf and icebergs. J. Phys. 
Oceanogr., 9, 189-198. 

Hellmer, H. H., and D. J. Olbers, 1989: A two-dimensional model 
for the thermohaline circulation under an ice shelf. Antarct. Sci., 
1, 325-336. 

, S. S. Jacobs, and A. Jenkins, 1998: Oceanic erosion of a floating 

Antarctic glacier in the Amundsen Sea. Ocean , Ice . and Atmo- 
sphere: Interactions at the Antarctic Continental Margin, S. S. 
Jacobs and R. F. Weiss, Eds., Antarctic Research Series, Vol. 
75, Amer. Geophys. Union, 83-99. 

Holland, D. M., and A. Jenkins, 1999: Modeling thermodynamic ice- 
ocean interactions at the base of an ice shelf. J. Phys. Oceanogr., 
29, 1787-1800. 

Huang, R. X., 1993: Real freshwater flux as a natural boundary con- 
dition for the salinity balance and thermohaline circulation 
forced by evaporation and precipitation. J. Phys. Oceanogr., 23, 
2428-2446. 

Jenkins, A., 1991: A one-dimensional model of ice shelf-ocean in- 
teraction. J. Geophys. Res., 96, 20 671-20 677. 

, and A. Bombosch, 1995: Modeling the effects of frazil ice 


crystals on the dynamics and thermodynamics of Ice Shelf Water 
plumes. J. Geophys. Res., 100, 6967-6981. 

Josberger, E. G., 1983: Sea ice melting in the marginal ice zone. J. 
Geophys. Res., 88, 2841-2844. 

Maykut, G. A., and M. G. McPhee, 1995: Solar heating of the Arctic 
mixed layer. J . Geophys. Res., 100, 24 691-24 703. 

McPhee, M. G., 1983: Turbulent heat and momentum transfer in the 
oceanic boundary layer under melting pack ice. J. Geophys. Res., 
88, 2827-2835. 

, 1990: Small-scale processes. Polar Oceanography, Part A: 

Physical Science, W. O. Smith Jr., Ed., Academic Press, 287- 
334. 

, G. A. Maykut, and J. H. Morison, 1987: Dynamics and ther- 
modynamics of the ice/upper ocean system in the marginal ice 
zone of the Greenland Sea. J. Geophys. Res., 92, 7017-7031. 

, C. Kottmeier, and J. H. Morison, 1999: Ocean heat flux in the 

central Weddell Sea during winter. J. Phys. Oceanogr., 29, 1 1 66— 
1179. 

Mellor, G. L., M. G. McPhee, and M. Steele, 1986: Ice-seawater 
turbulent boundary layer interaction with melting and freezing. 
J. Phys. Oceanogr., 16, 1829-1846. 

Millero, F. J., 1978: Annex 6: Freezing point of seawater. Eighth 
Report of the Joint Panel of Oceanographic Tables and Standards, 
UNESCO Tech. Paper Mar. Sci. 28, 29-31. 

Nicholls, K. W., and A. Jenkins, 1993: Temperature and salinity be- 
neath Ronne Ice Shelf, Antarctica. J. Geophys. Res., 98, 22 553- 
22 568. 

N0st, O. A., and A. Foldvik, 1994: A model of ice shelf-ocean in- 
teraction with application to the Filchner- Ronne and Ross Ice 
Shelves. J. Geophys. Res., 99, 14 243-14 254. 

Scheduikat, M., and D. J. Olbers, 1990: A one-dimensional mixed 
layer model beneath the Ross Ice Shelf with tidally induced 
vertical mixing. Antarct. Sci., 2, 29-42. 

Steele, M., G. L. Mellor, and M. G. McPhee, 1989: Role of the 
molecular sublayer in the melting or freezing of sea ice. J. Phys. 
Oceanogr., 19, 139-147. 


Oceanographic conditions beneath Ronne Ice Shelf: 
a comparison between model and field data 


Adrian Jenkins \ David M. Holland* and Keith W. Nicholls 

1 British Antarctic Survey, Natural Environment Research Council , 
High Cross , Madingley Road , Cambridge, CBS OET, U.K. 

2 Courant Institute of Mathematical Sciences , New York University, 
251 Mercer Street, New York City, NY 10012, U.S.A. 


Introduction 

Within the cavity beneath Filchner-Ronne Ice Shelf (FRIS), High Salinity Shelf Water 
(HSSW) is transformed into Ice Shelf Water (ISW) by a combination of melting and 
freezing at the ice shelf base and mixing within the water column. Evidence of these 
processes is abundant, both in the distribution of meteoric and marine ice in the ice shelf 
and in the properties of the water masses occupying the continental shelf. Quantitative 
information on the rate of water mass transformation is harder to come by, and more 
elusive still are indications of the rate-controlling factors. Is the circulation primarily 
“pulled” by the injection of meltwater deep beneath the ice shelf, or “pushed” by the 
production of dense waters at the ice front? Are there dynamical constraints on the rate 
of exchange at the ice front? Answers to these questions are highly relevant to the 
problem of how sensitive the sub-ice cavities are to externally forced climatic change. 

During the mid 1990’s a series of instrument strings were deployed through FRIS, in and 
near to the Ronne Depression (Nicholls and Makinson, 1998). Temperature and current 
records from the two sites within the depression showed a clear seasonality in the 
strength of the HSSW inflow; evidence of the importance of external forcing even deep 
within the cavity. The nature of the seasonal cycle, in particular the rapid increase in the 
HSSW flux associated with wintertime freezing north of the ice front, was further 
suggestive of some dynamical control on the inflow. 

In this paper we attempt to explain some of the key features of the seasonal signal using 
the results of a three-dimensional numerical model of ocean circulation within the FRIS 
cavity. We use a version of the Miami Isopycnic Coordinate Ocean Model (MICOM) 
that we have adapted for sub-ice-shelf domains (Holland and Jenkins, in press). The 
domain we use in this study is identical to that described by Jenkins and Holland (2000). 
The forcing is also almost identical, except that we have slightly modified the mid-winter 
surface salinity (Figure 1) to reflect the lower salinities that are currently generated over 
Berkner Bank (Nest and 0sterhus, 1998). The model was run for 10 years. 




Time series data from the cavity 

We focus on data from sites 2 and 3 (Figure 2). Two year, asynchronous temperature 
records from instruments suspended 30 to 80 m above the seabed at each site are shown 
in Figure 3. Both show relatively long periods with high temperatures, approaching the 
surface freezing point at site 2, followed by a gradual cooling and abrupt warming. 
Nicholls and Makinson (1998) interpreted the cooling trend as a weakening of the inflow 
to Ronne Depression, resulting from cessation of HSSW production at the ice front in 
summer. The warming was interpreted as the arrival of newer HSSW triggered by the 
onset of sea ice production and overturning of the water column north of the ice front. 



Figure 2: Sites mentioned in the text and bathymetry (m) throughout the model domain. 




Figure 4: Time series of model temperatures (°C) from grid points close to the seabed at 
the locations shown in Figure 2. 






The changes in temperature were assumed to reflect changes in the thickness of a lower 
layer of HSSW having relatively constant properties. The cooling was assumed to reflect 
the sinking of the thermocline past the level of the instrument. Of particular interest are 
the phasing and abruptness of the warming, associated with a thickening of the lower 
layer. The delay in arrival of the new HSSW was assumed to be related to the advection 
timescale from the ice front to the locations of the moorings, but the asymmetry between 
the abrupt warming and more gradual cooling was not explained. 

Figure 4 shows the temperature signal from the last two years of the model run at the grid 
points corresponding to the positions of the instrumental records described above. While 
some features of the model results, such as the amplitude of the variability at site 2 and 
the mean temperature at site 3, do not match the observations, other features, such as the 
asymmetry of the cycles and relative phasing of the warming at the two sites, show 
quantitative similarities with the instrumental records. The absolute timing of the model 
temperature minima is about 2 months out, but this is not surprising given the simple 
symmetrical forcing that is applied north of the ice front. Given that symmetry in the 
forcing, the asymmetry in the model response is a particularly interesting feature, the 
source of which may give insight into the processes affecting the observations. 



Figure 5: Time series of model salinity for grid points at the surface and near the seabed 
at the ice front location shown in Figure 2. 


We illustrate the impact of the surface forcing at the ice front with time series of model 
salinity from 0 and 500 m at a point near the centre of Ronne Depression (Figure 2). The 
temperature at both levels shows relatively little variability, being held close to the 
surface freezing point for most of the time. At the surface the salinity deviates only 
marginally from the sinusoidal restoring cycle, but at 500 m an asymmetry is already 
apparent (Figure 5). The salinity at any depth in the water column is only forced directly 
by the atmosphere when the surface mixed layer has deepened sufficiently. At 500 m 
direct atmospheric forcing is felt only between mid-July and the end of October. After 
this brief period of HSSW renewal, the springtime retreat of the mixed layer leaves a 
water column consisting of a stack of isopycnic layers. As the HSSW spreads away from 
its site of formation, the layers thin and, at the 500 m level, the salinity falls with the 
sinking of each layer interface past this depth. Convection once again reaches this level 
when the rising surface salinity of the following winter becomes equal to, and eventually 
exceeds, the salinity at depth. 

The impact of the relatively sharp onset of HSSW production is accentuated beneath the 
ice shelf by a coincident strengthening of the inflow to the cavity (Figure 6). While the 
flow is approximately parallel to the ice front during summer and autumn, as convection 
north of the ice front approaches full depth, the flow rapidly turns beneath the ice front. 
This is followed by a period of steady flow into the cavity, lasting until the mixed layer 
starts to retreat, at which time a gradual decrease in the inflow velocity begins. The result 
is the advection of a blob of new HSSW into the cavity, the front of which is particularly 
sharp. As the blob moves south along Ronne Depression, the domed isopycnals tend to 
slump, so that the signature of its arrival becomes progressively more diffuse. This 
thickening and thinning of the HSSW layer is the origin of the temperature variations 
shown in Figure 4. 



Figure 6: Time series of depth-mean velocity (cm s' 1 ) perpendicular to the ice front at the 
grid point shown in Figure 2. Negative values indicate flow into the cavity. 


Summary 


We have explored the propagation of the seasonal signature of HSSW production into the 
cavity beneath Filchner-Ronne Ice Shelf using MICOM. The focus of our study has been 
Ronne Depression, from where instrumental records are available for comparison. 
Agreement between model and observations is not perfect, but the phasing and shape of 
the annual cycle at two points along the depression is reasonably well simulated. Correct 
phasing implies that the speed of the inflowing HSSW is about right in the model. The 
asymmetric response of the model is intriguing, given the simple, symmetric forcing. We 
find that a combination of two factors is responsible for this. Production of HSSW only 
occurs for a relatively short period as the surface salinity approaches its maximum. At 
other times the dc:p parts of the depression are isolated from the surface buoyancy 
forcing. As HSSW production begins each year, there is a rapid transition from flow that 
is predominantly along the ice front, to flow that is directed almost perpendicular to the 
ice front. The new HSSW propagates into the cavity with a sharp front, that becomes 
gradually more diffuse with distance from the ice front. The seasonal forcing is felt 
throughout the model cavity and this, combined with the weak flow reported by Holland 
and Jenkins (in press) for an ideal cavity with no external forcing, suggests that our sub- 
ice version of MICOM is more “pushed” than “pulled”. 
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Abstract: Although it is isolated from direct atmospheric forcing, the ocean cavity 
beneath Filchner-Ronne Ice Shelf responds to the seasonal cycle of water mass 
production in the Weddell Sea. We have investigated the propagation of newly-formed 
shelf waters into the cavity using a version of the Miami Isopycnic Coordinate Ocean 
Model. The geographical focus of our study has been Ronne Depression, a seabed trough 
that crosses the Ronne Ice Front at its western end, and from where instrumental records 
are available for comparison with model output. Rapid warming is associated with the 
arrival of the new shelf waters at the instrumented sites beneath the ice shelf, and is 
followed by a more gradual cooling. The inflow to the cavity is regulated by the density 
structure in the water column to the north of the ice front, which in turn is controlled by 



the seasonal advance and retreat of the mixed layer. Within the Depression, the inflow 
peaks during the short period around mid-winter when convection reaches full-depth and 
the densest waters are generated. Once the surface density starts to decline, the deeper 
parts of the Depression are isolated from the surface buoyancy forcing. Dynamic 
adjustment of the restratified water column leads to a gradual fall in the bottom salinity 
and an associated weakening of the inflow. 

Introduction 

Important water mass conversions take place beneath Antarctica’s floating ice 
shelves. In the Ross and Weddell seas. High Salinity Shelf Water (HSSW) flows beneath 
the two largest ice shelves. There it is transformed into Ice Shelf Water (ISW), one of the 
precursors of Antarctic Bottom Water (AABW), by a combination of melting and 
freezing at the ice shelf base and mixing within the water column. HSSW forms through 
the enhancement of the shelf water salinity in regions where sea ice grows and is 
subsequently transported away, allowing further ice growth. While the net surface 
freshwater flux is determined by the wintertine atmospheric forcing, the properties and 
quantity of the resulting HSSW are also influenced by the rate at which the newly-formed 
shelf waters spread from their source region. Similarly, while the exchange of heat and 
freshwater at the ice shelf base effects the conversion of HSSW into ISW, the properties 
and quantity of ISW formed are determined by the rate at which HSSW ventilates the 
sub-ice cavity. Hence, the process of cavity ventilation has implications both for the 
properties of the abyssal ocean and for the mass budget of the Antarctic Ice Sheet, and 
represents one way in which the long-term evolution of ice sheet and ocean are linked. 



This paper discusses the details of the ventilation process, as revealed by observations 
and the results of a recently developed model of the circulation beneath Filchner-Ronne 
Ice Shelf (FRIS). 

Early models of sub-ice-shelf circulation [MacAyeal, 1984, 1985; Hellmer and 
Olbers , 1989; Jenkins , 1991] assumed a free exchange of waters across the ice front, 
driven by the density difference between the inflow and outflow. Later simulations 
[Hellmer and Olbers, 1991; Hellmer and Jacobs, 1992] illustrated how the ice shelf and 
seabed topography could combine to leave regions of the cavity relatively isolated from 
the waters beyond the ice front. The first three-dimensional ocean general circulation 
models to be applied to sub-ice-shelf domains [Determann and Gerdes, 1994; Grosfeld et 
al., 1997] further highlighted the role of topography in guiding the circulation. Strong 
flows occurred parallel to water column isopachs, with little cross-isopach transport. In 
these models, the step in water column thickness at the ice front represented a 
considerable barrier to flow, such that only diffusive exchange occurred between the 
cavity and open ocean, unless the topography of the seabed and ice shelf base combined 
to allow continuity of isopachs across the ice front. The implication of these results is 
that the cavity circulation may be relatively insensitive to the atmospheric forcing north 
of the ice front, and may therefore be insulated from climatic change. 

Observational data are sparse, particularly within the cavities, but suggest that 
outflows from beneath the Ross and Filchner-Ronne ice shelves each total ~1 Sv [Foldvik 
et al., 1985; Jacobs et al, 1992], At this rate the cavity waters could be completely 
renewed on a timescale of 5-10 years, similar to the residence times inferred from tracer 
studies of the emerging waters [Trumbore et al, 1991; Gammelsr0d et al., 1994]. During 



the mid 1990’s a series of instrument strings were deployed through FRIS in Ronne 
Depression, a seabed trough parallel to the western coast of the ice shelf (Figure 1). 
Temperature and current records from beneath the ice shelf showed a clear seasonality in 
the strength of the HSSW inflow [Nicholls and Makinson , 1998], providing evidence of 
the importance of external forcing even deep within the cavity. The nature of the 
seasonal cycle suggested that the dynamical control on the inflow was more complex 
than could be explained simply by the arrangement of water column isopachs. Current 
meters deployed at Site FR6, on the western flank of Ronne Depression just north of the 
ice front (Figure 1), confirmed this as a location of inflowing HSSW [Woodgate et al , 
1998]. A pronounced seasonal variation was found in the velocity perpendicular to the 
ice front. In contrast to the results of the three-dimensional models described above, 
these observations imply that much of the cavity is impacted by the seasonal atmospheric 
forcing north of the ice front, and so is likely to react relatively rapidly to changes in that 
forcing. 

Time series observations from Ronne Depression 

Two year, asynchronous temperature records from instruments suspended 30 to 
80 m above the seabed at sites 2 and 3 (Figure 1) are shown in Figure 2a, b. Both show 
relatively long periods with high temperatures, approaching the surface freezing point at 
Site 2, followed by a gradual cooling and abrupt warming. Nicholls and Makinson 
[1998] interpreted these changes as the propagation of a seasonal signal, associated with 
the wintertime production of HSSW at the ice front, along the Ronne Depression. 
Following their arguments the cooling trend at each site is a result of the weakening of 



the HSSW inflow to the cavity during the summer, while the abrupt warming is triggered 
by the onset of sea ice production north of the ice front and the inflow of new HSSW. 

Nicholls and Makinson [1998] assumed that the changes in temperature reflected 
changes in the thickness of a lower layer of HSSW that had relatively constant properties. 
Warming and cooling trends recorded by each instrument were assumed to indicate the 
rising and falling of a thermocline past the level of the sensor, while periods of relatively 
constant temperature were thought to indicate the complete immersion of the instrument 
in the HSSW layer. 

Of particular interest are the phasing and abruptness of the warming, associated 
with the thickening of the lower layer. The delay in arrival of the new HSSW was 
assumed to be related to the advection timescale from the ice front to the locations of the 
moorings. This gave an estimate of the inflow velocity that was consistent with estimates 
derived by other means. The asymmetry between the abrupt warming and more gradual 
cooling was not explained. 

This simple picture of the seasonal intrusion of newly-formed HSSW into the 
cavity is drawn into question by the observations at Site FR6, where currents were 
measured at three levels in the water column over a period of about 2 yr. The records 
indicate a flow that had a relatively constant magnitude throughout the year, but a 
seasonally varying direction. The component of the velocity perpendicular to the ice 
front (Figure 2c) peaked in mid-summer and fell close to zero by mid-winter. So the 
onset of surface freezing in the coastal polynya is associated with a weakening flow into 
the cavity at this point, and the inflow of new HSSW does not begin to build until late 
winter, after its apparent arrival at Site 2. 



The aim of this paper is to use a three-dimensional numerical model of ocean 
circulation within and to the north of the FRIS cavity to resolve this apparent 
discrepancy, and to help explain the character of the signals recorded beneath the ice 
shelf. 

The model 

We use a version of the Miami Isopycnic Coordinate Ocean Model (MICOM) that 
has been adapted for sub-ice-shelf domains [Holland and Jenkins , in press]. In this study 
our model grid covers the entire continental shelf of the southern Weddell Sea (Figure 1) 
at a resolution of 0.75° of longitude. We use 10 isopycnic layers with equal density steps 
between them, the density range being chosen to be consistent with that observed within 
the region covered by the model domain [Olbers et al. t 1992]. In addition there is a 
surface mixed layer, which has a spatially and temporally varying density. 

The model is forced by restoring the mixed layer temperature and salinity in the 
open ocean to a prescribed pattern. The mid-winter surface salinity field is defined to 
increase both to the south and to the west and to decrease with increasing water depth, 
giving a pattern that is broadly similar to that of the bottom salinity observed over the 
shelf during summer [Foldvik et ai , 1985]. In this way we ensure that the distribution 
and properties of the HSSW that we generate in the model approximate those found in 
nature. The surface temperature is defined to be equal to the salinity-dependent freezing 
point at all times, consistent with the presence of a perennial ice cover over much of the 
domain. The mid-summer surface salinity field has a similar pattern to that of mid- 
winter, and the variation between summer and winter extremes is a simple sinusoid, the 



amplitude of which increases to the south and west. This reflects the distribution of 
winter polynyas, which occur all along the Filchner-Ronne Ice Front, but are most 
extensive in the west [Markus et al., 1998], where persistent barrier winds drive the ice 
north. In the model there is no wind forcing anywhere, while beneath the ice shelf the 
only direct forcing is the melting and freezing that is computed at the ice-ocean interface. 
Initial conditions are set according to Olbers et al. [1992], with simple extrapolation to 
fill the region beneath the ice shelf. The northern and eastern boundaries of the domain 
are solid walls, where we restore to initial conditions. The model is run for 10 years 

Seasonal response of the cavity interior 

Figure 3a shows the temperature signal from the last two years of the model run at 
the grid points and depths corresponding to the positions of the instrumental records 
shown in Figure 2a,b. While some features of the model response, such as the amplitude 
of the variability at Site 2 and the mean temperature at Site 3, do not match the 
observations, other features, such as the asymmetry of the cycles and the relative phasing 
of the temperature minima at the two sites, show quantitative similarities with the 
instrumental records. The times of the model temperature minima are about 2 months 
late, but this is not surprising given the simple symmetrical forcing that is applied north 
of the ice front. Given that symmetry in the forcing, the asymmetry in the model 
response is a particularly interesting feature, the source of which we address in the next 
section. 

The process responsible for the variability shown in Figure 3 is precisely that 
envisaged by Nicholls and Makinson [1998]. HSSW flows along Ronne Depression as a 


dense current that hugs the eastern flank (Figure 4). It is thus guided past the locations of 
the two moorings and the seasonal renewal of the HSSW causes the dense layer to 
abruptly thicken. This signature becomes progressively delayed and attenuated as it 
propagates along the Depression. At the location of site 2 the arrival of the new HSSW is 
accompanied by a 3-4 fold increase in basal melting (Figure 3b). While there is also a 
distinctive seasonal variation in melting at site 3, the peak is not coincident with the 
arrival of the HSSW. The warming there seems to be too deep in the water column to 
have a direct impact on the ice shelf. The phasing of the melt signature is likely to be the 
result of a further delay while the wanning in the deep layers is advected to a site that is 
more favourable for upwelling, such as a grounding line, and may reflect a regional 
response to the increased heat supply to the mixed layer at the grounding line. 

Forcing and response at the ice front 

The net surface freshwater flux required to drive the mixed layer in Ronne 
Depression through its prescribed salinity cycle is illustrated in Figure 5a. This quantity 
is directly comparable to the growth rate of sea ice, estimates of which have been made 
for the entire coastal polynya of the southern Weddell Sea by Renfrew et al [in press]. 
For the period 1992 to 1998, they estimate that the freezing season started on average in 
mid- to late February and ended on average in early November. Compared with these 
results, our model freezing season is of approximately the correct length, but occurs 
between 1 and 2 months late. The lag in the freezing cycle presumably explains why the 
temperature variations at sites 2 and 3 in our model appear to lag the observations by a 
couple of months. Overall the forcing is rather weak, with an average freshwater flux of 



7.5 m yr' 1 over the freezing season compared with over 20 m yr' 1 estimated by Renfrew et 
al. [in press]. This may reflect the lack of wind forcing in the model, which would tend 
to drive HSSW northwards away from its site of formation in the coastal polynya. With 
the winds included in the forcing, we might anticipate that a larger freshwater flux would 
be required to restore the mixed layer to the prescribed salinity pattern. 

We illustrate the impact of the surface forcing at the ice front with time series of 
model salinity (Figure 5b) from 0 and 500 m depth at the grid point corresponding to Site 
FR6 (Figure 1). The temperature at both levels shows relatively little variability, being 
held close to the surface freezing point for most of the time. At the surface the salinity 
deviates only marginally from the sinusoidal restoring cycle, but at 500 m an asymmetry 
is apparent. The salinity at any depth in the water column is only forced directly by the 
atmosphere when the surface mixed layer has deepened sufficiently. At 500 m direct 
atmospheric forcing is felt only between mid-July and the end of October. After this 
brief period of HSSW renewal, the springtime retreat of the mixed layer leaves a water 
column containing significant horizontal and vertical density gradients. Isopycnal 
surfaces are initially domed around the regions where the densest HSSW has formed. As 
the HSSW spreads away from its site of formation the isopycnals slump and, at the 500 m 
level in Ronne Depression, the salinity falls with the sinking of each layer interface past 
this depth. Convection once again reaches this level when the rising surface salinity of 
the following winter becomes equal to, and eventually exceeds, the salinity at depth. 

The impact of the relatively sharp onset of HSSW production is accentuated 
beneath the ice shelf by a coincident strengthening of the inflow to the cavity. During the 
summer and autumn a steady current of -10 cm s 1 carries fresher water from east to west 



along the ice front, before turning north along the western coast (Figure 6a). At this time 
of year the current is confined to the surface layers. Geostrophic adjustment of the 
deeper layers leads to anti-cyclonic circulation around the regions of maximum salinity 
(Figure 6b). In the far west, a southerly flow carries HSSW into the cavity. This flow 
persists until the waters have been drained from the relatively shallow region to the north 
of FR6 (Figure 1). As convection deepens the mixed layer, more waters are drawn into 
the current that parallels the ice front until it extends over the full depth of the water 
column. This combined with the rising salinity in the shallower waters to the north of 
FR6, means that this current can no longer turn north along the coast, and is forced south 
beneath the ice shelf (Figure 6c). Melting at the ice front immediately restratifies the 
incoming water column, and almost unmodified HSSW drains south-west along Ronne 
Depression (Figure 6d). The melting contributes towards an outflow of ISW around 
58°W, which is present in summer also, but is more diffuse then. This outflow is a 
feature that is common to all observations of oceanographic conditions at the ice front 
[e.g. Foldvik at al., 1985; Gammelsr0d et al., 1994]. 

Averaged across the width of Ronne Depression the seasonal cycle in the inflow 
shows a dramatic saw-tooth shape (Figure 7a). The rapid increase in velocity 
perpendicular to the ice front begins as the surface forcing reaches its maximum and 
convection reaches full depth at FR6 (Figure 5). The warming at the location of Site 2 
follows about one month later (Figure 3). As the surface forcing decays, and the water 
column restratifies, the decrease in the inflow velocity is more gradual. There is a 
secondary peak in late summer, corresponding to the final stages of adjustment to the 
previous winter’s convection (Figure 6b). At the westernmost grid point along the ice 


front, the seasonal cycle in the inflow (Figure 7b) is dominated by this secondary summer 
peak, the timing and magnitude of which, at this point, agree favourably with those of the 
maximum inflow recorded at Site FR6. However, the mid-winter inflow clearly has the 
greatest impact on the cavity, and the sudden onset and more gradual cessation of this 
flow is the origin of the asymmetry seen in Figure 3a. 

The absence of any significant time lag between the observed and modelled peak 
inflow at the western end of the ice front probably results from the simplicity of the 
surface forcing (Figure 5a). Compared with observation [Renfrew et al . , in press] the 
freezing peak is too narrow. The longer, more gradual decline in the model freezing 
cycle means that the deep layers within the Depression begin their post-convection 
adjustment too soon. A premature weakening of the mid-winter inflow could also 
explain why the model temperature records (Figure 2b) fail to reproduce the period of 
sustained high temperatures that follow the abrupt warming seen in the observations 
(Figure 2a, b). 

Summary and Conclusions 

We have explored the propagation of the seasonal signature of HSSW production 
into the cavity beneath Filchner-Ronne Ice Shelf using MICOM. Agreement between 
model and observations is not perfect, but the phasing and shape of the annual 
temperature cycle at two points along Ronne Depression is reasonably well simulated. 
Correct phasing implies that the speed of the inflowing HSSW is about right in the 
model. The only direct observations of inflow to the cavity in this region come from the 
west of the Depression. Once again some of the details of the model results do not match 



the observations, but the timing and magnitude of the maximum velocity are simulated 
reasonably well. 

The response of the model cavity is governed by the density structure within the 
water column to the north of the ice front. Density reaches a maximum for a short period 
around mid- winter. This is when the most saline shelf waters are renewed and, 
simultaneously, flow into the cavity increases. These processes are forced directly by the 
surface freshwater flux, so their onsets are relatively rapid. The new HSSW propagates 
into the cavity with a sharp front that becomes gradually more diffuse with distance from 
the ice front. The inflow switches off more gradually as the water column north of the 
ice front is restratified, leaving the deeper layers once more isolated from the surface 
buoyancy forcing. This phase of decreasing salinity and generally weakening flow is 
governed by the dynamic adjustment of the deeper parts of the water column to the 
earlier forcing. Anti-cyclonic circulation around domed isopycnals tends to carry HSSW 
northwards in the middle of the Depression, but in the far west drainage into the cavity 
continues and even increases in strength into the summer. The impact of the seasonal 
forcing is seen throughout the model cavity, although the phasing is different in different 
geographical locations and at different levels in the water column. 

The above analysis would suggest that the inflow recorded at Site FR6 is of 
secondary importance for the region further south and west along Ronne Depression, 
when compared with inflows that occur immediately to the east of the mooring site at 
other times of the year. This might seem a bold statement, if it were based solely on the 
results from a crudely-forced, relatively low resolution model, with smoothed and 
possibly erroneous seabed topography. However, we should re-emphasise that the 



observations at sites 2 and 3 are difficult to explain if the main inflow of HSSW to the 
Ronne Depression peaks in mid-summer. We contend that the model demonstrates a 
plausible mechanism through which the records from the cavity and that at FR6 may be 
reconciled. 

The exact processes by which the signature of seasonal forcing north of the ice 
front is transmitted to the interior of the model cavity will be an abstraction of those 
processes that operate in nature. However, the realistic features of the model response in 
Ronne Depression lends us confidence in drawing some important, general conclusions. 
It is clear that even remote parts of the sub-ice-shelf cavity are impacted by external 
forcing on sub-annual time-scales. This is because the exchange of waters across the ice 
front is controlled not just by the arrangement of water column isopachs, but by the three- 
dimensional distribution of density within the water column. We should therefore also 
anticipate a cavity response to inter-annual variations in sea ice production and to longer- 
term climatic trends. Such a response has implications for the mass balance of the ice 
shelf and the production of AABW. 

Acknowledgements: This work has received financial support from the U.S. National 
Aeronautics and Space Administration (grant NAG-5-4028). 

References 

Determann, J., and R. Gerdes, Melting and freezing beneath ice shelves. Implications 
from a three-dimensional ocean-circulation model, Ann . GlacioL, 20 , 413—419, 



Foldvik, A., T. GammeIsr0d, and T. T0rresen, Circulation and water masses on the 
southern Weddell Sea shelf, in Oceanology of the Antarctic Continental Shelf, 
Antarct. Res. Ser., vol. 43, edited by S. S. Jacobs, pp. 5-20, AGU, Washington, 
D.C., 1985. 

Gammelsrpd, T„ A. Foldvik, O. A. Npst, 0. Skagseth, L. G. Anderson, E. Fogelquist, K. 
Olsson, T. Tanhua, E. P. Jones, and S. 0sterhus, Distribution of water masses on 
the continental shelf in the Southern Weddell Sea, in The Polar Oceans and Their 
Role in Shaping the Global Environment, edited by O. M. Johannesen, R. D. 
Muench, and J. E. Overland, pp. 159-176, AGU, Washington, D.C., 1994. 

Grosfeld, K., R. Gerdes, and J. Determann, Thermohaline circulation and interaction 
between ice shelf cavities and the adjacent open ocean, J. Geophys. Res., 102, 
15,595-15,610, 1997. 

Hellmer, H. H., and S. S. Jacobs, Ocean interactions with the base of the Amery Ice 
Shelf, Antarctica, J. Geophys. Res., 97, 20,305-30,317, 1992. 

Hellmer, H. H., and D. J. Olbers, A two-dimensional model for the thermohaline 
circulation under an ice shelf, Antarct. Sci., I, 325-336, 1989. 

Hellmer, H. H., and D. J. Olbers, On the thermohaline circulation beneath the Filchner- 
Ronne Ice Shelves, Antarct. Sci., 3, 433—4-42, 1991. 

Holland, D. M., and A. Jenkins, Adaptation of an isopycnic coordinate ocean model for 
the study of circulation beneath ice shelves, Mon. Wea. Rev., in press. 

Jacobs, S. S., H. H. Hellmer, C. S. M. Doake, A. Jenkins, and R. M. Frolich, Melting of 
ice shelves and the mass balance of Antarctica, J. Glaciol., 130, 375-387, 1992. 


Jenkins, A., A one-dimensional model of ice shelf-ocean interaction, J. Geophys. Res., 
96, 20,671-20,677, 1991. 

MacAyeal, D. R., Thermohaline circulation below the Ross Ice Shelf: A consequence of 
tidally induced vertical mixing and basal melting, J. Geophys. Res., 89, 597-606, 
1984. 

MacAyeal, D. R., Evolution of tidally triggered meltwater plumes below ice shelves, in 
Oceanology of the Antarctic Continental Shelf, Antarct. Res. Ser., vol. 43, edited 
by S. S. Jacobs, pp. 133-143, AGU, Washington, D.C., 1985. 

Markus, T., C. Kottmeier, and E. Fahrbach, Ice formation in coastal polynyas in the 
Weddell Sea and their impact on oceanic salinity, in Antarctic Sea Ice: Physical 
Processes, Interactions and Variability, Antarct. Res. Ser., vol. 74, edited by, M. 
O. Jeffries, pp. 273-292, AGU, Washington, D.C., 1998. 

Nicholls, K. W., and K. Makinson, Ocean circulation beneath the western Ronne Ice 
Shelf, as derived from in situ measurements of water currents and properties, in 
Ocean, Ice, and Atmosphere: Interactions at the Antarctic Continental Margin, 
Antarct. Res. Ser., vol. 75, edited by S. S. Jacobs and R. F. Weiss, pp. 301—318, 
AGU, Washington, D.C., 1998. 

Nost, O. A. and S. Osterhus, Impact of grounded icebergs on the hydrographic conditions 
near the Filchner Ice Shelf, in Ocean, Ice, and Atmosphere: Interactions at the 
Antarctic Continental Margin, Antarct. Res. Ser., vol. 75, edited by S. S. Jacobs 
and R. F. Weiss, pp. 267-284, AGU, Washington, D.C., 1998. 

Olbers, D. J., V. Gouretski, G. SeiB, and H. Schroter, Hydrographic Atlas of the Southern 
Ocean, Alfred- Wegener-Institute, Bremerhaven, Germany, 1992. 



Renfrew, I. A., J. C. King, and T. Markus, Coastal polynyas in the southern Weddell Sea: 
Variability of the surface energy budget, J. Geophys. Res., in press. 

Trumbore, S., S. S. Jacobs, and W. Smethie, Jr., Chloroflourocarbon evidence for rapid 
ventilation of the Ross Sea, Deep-Sea Res., 38, 845-870, 1991. 

Woodgate, R. A., M. Schroder, and S. 0sterhus, Moorings from the Filchner Trough and 
the Ronne Ice Shelf Front: Preliminary results, in Filchner-Ronne Ice Shelf 
Programme, Report No. 12, edited by H. Oerter, pp. 85-90, Alfred-Wegener- 
Institute for Polar and Marine Research, Bremerhaven, 1998. 



V ! 


70S 

71S 

72S 

73S 

74S 

75S 

76S 

77S 

78S 

79S 

80S 

81S 

82S 



80W 75W 70W 65W 60W 55W SOW 45W 40W 35W 30W 25W 20W 


Figure 1: Model domain covering the continental shelf of the southern Weddell Sea. 
Contour lines indicate seabed depth (m). The region to the south and west of the bold 
line is covered by Filchner-Ronne Ice Shelf. Labelled solid circles indicate grid points 
that correspond to the location of moored instruments, the records from which are 
discussed in the text. The dashed line passing through Site 2 indicates the section plotted 
in Figure 4. The box shown around Site FR6 indicates the area used to characterise the 
surface forcing over Ronne Depression (Figure 5), the seabed trough in which all the 
moorings are located. 
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Figure 2: Time series of potential temperature recorded near the seabed at a) Site 2, and 
b) Site 3. c) Time series of currents perpendicular to the ice front recorded at Site FR6, 
Negative values indicate flow into the cavity. 
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Figure 3: a) Time series of model potential temperature from two of the locations shown 
in Figure 1. The solid line is an average over the depth range 750 to 850 m below sea 
level from Site 2, while the dashed line is an average over 1 100 to 1200 m at Site 3. b) 
Time series of modelled melt rate at the ice shelf base for Site 2 (solid line) and Site 3 




(dashed line). 





Figure 4: Cross-sections of model potential temperature along the line shown in Figure 1 : 
a) in late July, before the arrival of the HSSW inflow, and b) in late August after the 
arrival of the inflow. The vertical line indicates the location of Site 2. Light shading 
denotes the ice shelf, dark shading the seabed. 



FEB APR JUN AUG OCT DEC FEB APR JUN AUG OCT DEC 


f 

# 34.7 



FEB APR JUN AUG OCT DEC FEB APR JUN AUG OCT DEC 

Month of year 


Figure 5: a) Time series of net surface freshwater flux, averaged over the box shown by 
the dotted outline in Figure 1. Positive values indicate removal of freshwater from the 
ocean, the model analogue of surface freezing, b) Time series of model salinity for the 
grid point labelled Site FR6 in Figure 1 . The solid line indicates surface values, while the 
dashed line indicates values at a depth of 500 m. 
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Figure 6: Current vectors showing circulation in a) the mixed layer in late January, b) 
layer 9 in late January, c) the mixed layer in late September and d) layer 9 in late 
September. The bold line indicates the ice front, while light shading in b) and d) shows 
where layer 9 does not exist. Scale (cm s' 1 ) is indicated by the arrows at the bottom right 


of each panel. 






THE WEST ANTARCTIC ICE SHEET: BEHAVIOR AND ENVIRONMENT 
ANTARCTIC RESEARCH SERIES, VOLUME 77, PAGES 237-256 


A REVIEW OF PINE ISLAND GLACIER, WEST ANTARCTICA: 
HYPOTHESES OF INSTABILITY VS. OBSERVATIONS OF CHANGE 

David G. Vaughan 1 , Andrew M. Smith 1 , Hugh F. J. Corr', Adrian Jenkins 1 , Charles R. Bentley 2 , 
Mark D. Stenoien 2 , Stanley S. Jacobs 3 , Thomas B. Kellogg 4 , Eric Rignot 5 and Baerbel K. Lucchitta 






The Pine Island Glacier ice-drainage basin has been cited as the part of the 
West Antarctic ice sheet most prone to substantial retreat on human time-scales. 
Here we review the literature and present new analyses showing that this ice- 
drainage basin is glaciologically unusual. Due to high precipitation rates near the 
coast, Pine Island Glacier basin has the second highest balance flux of any extant 
ice stream or glacier. Well-defined tributaries flow at intermediate velocities 
through the interior of the basin and have no regions of rapid velocity increase. The 
tributaries coalesce to form Pine Island Glacier which has characteristics of outlet 
glaciers (e.g. high driving stress) and of ice streams (e.g. shear margins bordering 
slow-moving ice). The glacier flows across a complex grounding zone into an ice 
shelf. There, it comes into contact with warm Circumpolar Deep Water which fuels 
the highest basal me It- rates yet measured beneath an ice shelf. The ice front 
position may have retreated within the past few millennia but during the last few 
decades it appears to have shifted around a mean position. Mass balance 
calculations of the ice-drainage basin as a whole show that there is currently no 
measurable imbalance, although there is evidence that some specific areas within 
the basin are significantly out of balance. The grounding line has been shown to 
have retreated in recent years. The Pine Island Glacier basin is clearly important in 
the context of the future evolution of the West Antarctic ice sheet because 
theoretically, it has a high potential for change and because observations already 
show change occurring. There is, however, no clear evidence to indicate sustained 
retreat or collapse over the last few decades. 


1. INTRODUCTION 

The West Antarctic ice sheet (Figure 1) drains into the 
Southern Ocean by three main routes; through the ice 
streams on the Siple Coast into the Ross Ice Shelf, through 
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the glaciers and ice streams feeding Ronne Ice Shelf, and 
through the glaciers which debouch, either directly or 
through small ice shelves, into the Bellingshausen and 
Amundsen seas. While the dynamics of the Siple Coast ice 
streams have been studied under the West Antarctic Ice 
Sheet Initiative (WAIS), and those feeding Ronne Ice 
Shelf have been studied under the auspices of the Filchner- 
Ronne Ice Shelf Programme (FRISP), there has been no 
coordinated effort to understand the dynamics of glaciers 
feeding the Bellingshausen and Amundsen seas. 
Consequently, this area is seldom visited and its glaciology 
is poorly understood. 

The largest glaciers in this sector are Pine Island 
Glacier and Thwaites Glacier. Both transport ice from the 
interior of the West Antarctic ice sheet to the Amundsen 
Sea. In terms of the mass of snow accumulating in their 
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Fig. 1. Location map. Shaded region shows area covered by 
Figures 2-7. 

catchment basins* Pine Island and Thwaites glaciers are 
respectively, the second and fifth most active basins in 
Antarctica [Vaughan and Bamber, 1998]. Pine Island 
Glacier alone accounts for around 4% of the outflow from 
the entire Antarctic Ice Sheet. The ice-drainage basins that 
feed these glaciers rest on beds as much as 2500 m below 
sea level, perhaps the deepest in Antarctica, and some 
authors have suggested that this in itself implies a great 
potential for rapid collapse [e.g. Fastook, 1984; Thomas , 
1984]. 

Together, Pine Island and Thwaites glaciers may be key 
to the future evolution of the West Antarctic ice sheet, but 
in this review, we concentrate on the Pine Island Glacier 
basin alone. We do this because, in addition to theories of 
instability, there is a growing body of observations of 
change and unsteady flow there. After some introductory 
notes we consider each of the component parts of the basin 
in turn. We then consider the interactions between the 
basin and sea into which it flows. We assess the evidence 
for both long-term and recent changes in the ice cover of 
the region. Finally, we consider how relevant those 
observations may be to models which have predicted that 
Pine Island Glacier might be particularly prone to collapse. 

1 . 1 Introductory notes on nomenclature 

1.1.1 Glaciers , ice streams and outlet glaciers . The 
terms glacier , ice stream and outlet glacier are often 


loosely applied in the scientific literature as well as the 
non-specialist press. Here we use widely accepted 
definitions; ice streams , being areas of fast-moving ice 
sheet bounded by slower moving ice; outlet glaciers, being 
fast-moving ice bounded by nunataks or mountain ranges 
[ Bentley , 1987; Swithinbank, 1954]; and glaciers , being a 
generic term for a “mass of snow and ice continuously 
moving from higher to lower ground, or if afloat, 
spreading continuously” [Armstrong et al, 1973]. The 
distinction between ice streams and outlet glaciers 
“becomes rather hazy in practice” [Bentley, 1987], and is 
particularly acute in this case, as Pine Island and Thwaites 
glaciers share most of the dynamical characteristics of 
pure ice streams and need not be considered as inherently 
different. 

1.1.2 FLating portion of Pine Island Glacier. Hughes 
[ 1 980] stated that neither Pine Island Glacier nor Thwaites 
Glacier are “buttressed by a confined and pinned ice 
shelf’. Stuiver et al [1981] also stated that they “are 
unimpeded by an ice shelf’. At that time, the position of 
the grounding line of Pine Island Glacier was poorly 
mapped, and Pine Island Glacier was assumed to calve 
directly into Pine Island Bay. Airborne radar sounding 
soon revealed that the seaward -80 km of the glacier was 
indeed floating [Crabtree and Doake, 1982]. The original 
idea that these glaciers are dynamically different from 
others in West Antarctica has, however, persisted. 
Subsequent authors support the original notion, that Pine 
Island Glacier does not debouch through “an ice shelf’ 
[Kellogg and Kellogg, 1987], “a substantial ice shelf’ 
[Jenkins etal, 1997], or “a large ice shelf’ [ Rignot , 1998]. 
Taking the widely accepted definition of an ice shelf, a 
“floating ice sheet of considerable thickness attached to a 
coast " [Armstrong et al , 1973], it is clear that the floating 
portion of Pine Island Glacier, together with the adjacent 
floating areas, do constitute an ice shelf. Whether the ice 
shelf is a significant dynamic control on the glacier is, 
however, still an open question. 

1.1.3 Pine Island Bay . Strictly, Pine Island Bay is the 
bay (approximately 75 x 55 km) into which flows the ice 
from Pine Island Glacier [Alberts, 1981], although the term 
is sometimes applied to a rather wider area. 

1.1.4 West Antarctic ice sheet. The West Antarctic ice 
sheet is not an officially-recognised placename. We follow 
the accepted usage, meaning the term to refer to the ice 
sheet that covers West Antarctica, but excluding the 
Antarctic Peninsula. 

1.2 Introductory note on meteorology 

Since the 1960s, it has been widely recognised that the 
coastal portions of West Antarctica bordering on the 
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Fig. 2. Map of the area around the ice-drainage basin of Pine Island Glacier showing sites of 
measured surface mass balance (black dots), and interpreted grid of surface mass balance derived 
from field measurements and passive microwave satellite data from Vaughan et al [1999]. 




Amundsen and Bellingshausen seas experience high 
precipitation rates compared with the rest of Antarctica 
[Shimizu, 1964]. These rates are matched only on the 
Antarctic Peninsula and around the coast of Wilkes Land 
[Giovinetto, 1964]. A recent compilation of net surface 
mass balance derived from in situ measurements and 
satellite data [Vaughan et al., 1999] is shown in Figure 2. 
It agrees broadly with earlier estimates [eg. Giovinetto and 
Bentley , 1985], but shows an increased level of detail. 

Meteorologically, the high precipitation rate in this 
sector results from synoptic-scale cyclones which travel 
around the Antarctic in the circumpolar trough. The trough 
is deepest over the Amundsen Sea, and many synoptic- 


scale cyclones come ashore here, producing considerable 
precipitation in the coastal zone. The effect is seen in 
moisture transport calculations [Bromwich, 1988] and in 
precipitation fields derived from General Circulation 
Models [e g., Connolley and KingJ996\. Precipitation may 
be particularly high during winter months when the 
circumpolar trough moves south and more cyclones track 
across the coastal region [Jones and Simmonds , 1993], 
There are no meteorological stations in the 
Bellingshausen/Amundsen Sea sector of West Antarctica. 
Consequently, direct measurements of decadal climate 
change (or stasis) have yet to be reported from either the 
Pine Island Glacier or Thwaites Glacier basins, although it 
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Fig. 3. Map of the area around Pine Island Glacier showing surface-elevation contours (meters 
above sea level) derived from Geodetic Mission of the ERS-1 satellite on a 5-km grid [Bamber 
and Bindschadier, 1997] Ice shelves are shaded 


is possible that surface elevation changes (see Section 2.7) 
do reflect recent anomalous precipitation rates [Wingham 
etal. , 1998]. Offshore, a reduction in sea-ice extent in both 
the Amundsen and Bellingshausen seas has been noted in 
all seasons over the two decades prior to 1 995 [ Jacobs and 
Comiso, 1997]. This perhaps reflects a change in surface 
temperatures. 

2. THE INTERIOR ICE-DRAINAGE BASIN 

The interior of ice-drainage basins are sometimes 
viewed as cisterns, which passively accumulate ice and 
then supply it to the glacier (or ice stream) at whatever rate 
the glacier can transport it away. In this section, however, 
we present evidence that flow in the Pine Island Glacier 
basin is far from homogeneous. There is no clear 
distinction between ice-sheet and glacier flow and flow in 


the basin may have a strong influence on the overall 
configuration and glacier activity. 

2.1 Delineation of the ice-drainage basin 

Field-measurements of surface elevation in the Pine 
Island Glacier basin are few, but altimetry from the ERS- 1 
satellite is available to 81.9°S, which includes the entire 
basin. This altimetry has been used to create several high- 
resolution Digital Elevation Models (DEMs) of Antarctica 
[eg., Bamber and Bindschadier, 1997; Legresyand Remy, 
1997 , Stenoien, 1998; Liu etal, 1999], 

Here we use a 5 km-resolution ERS- 1 -derived DEM 
[Bamber and Bindschadler, 1997] (Figure 3) to delineate 
the Pine Island Glacier basin and its neighbours. For 
comparison, a 200-m resolution DEM [Liu et al. , 1 999] 
was also used to delineate the Pine Island Glacier basin 
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Fig 4. Map of ice-drainage basins in the vicinity of Pine Island Glacier. Solid lines show basm 
boundaries derived from Bomber and B.ndschodler [1997], Dotted line shows the Pme Island 
Glacier basin derived from Liu e, al. [1999], Major glaciers are numbered and their d.rect.on of 
flow indicated by arrows; I. Pine Island Glacier, 2. Thwa.tes Glacier, 3. Evans Ice Stream, A 
Carlson Inlet, 5. Rutford Ice Stream, 6. Institute Ice Stream. Rock outcrops are hatched with the 
major mountain ranges numbered (7. Hudson Mountains, 8^ Jones Mountains, 9. Ellsworth 
Mountains) I and II indicate regions for which Stenoien [ 1 998 j presented septate mass-ba ance 
calculations. The coastline is derived from the Antarctic Digital Database [SCA , 1 

shelves are shaded. 


alone (Figure 4). The method used to produce the DEM is 
described in detail by Vaughan etal. [1999]. In summary, 
we identified segments of grounding line and then 
delineated the basins that feed them by tracing the line of 
steepest ascent inland as far as the ice divide. This 
procedure was limited to the grounded ice sheet as it 
assumes that ice-flow is parallel to the surface slope. 

Table 1 shows the area of the Pine Island Glacier 
drainage basin measured from the above delineation, using 
an equal area projection, together with earlier estimates for 
comparison. There was considerable disagreement between 
the early estimates. While the more recent ones that rely on 


ERS-1 data have reduced the uncertainty, some still 
remains, presumably resulting from the different methods 
of analysis. For this review, we use an average of the three 
most recent values (165,000±7000 km 3 ) as a reasonable 
estimate of the basin area, although we accept that this may 
be improved in future by further analysis. 

2.2 Shape of the catchment basin 

The shape of the Pine Island Glacier basin (Figure 4) is 
similar to earlier delineations. It consists of two lobes, one 
immediately upstream of Pine Island Glacier and another. 
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Table 1 . Estimates of ice-drainage basin area and balance flux for the Pine Island Glacier basin. 


Basin Area 

Balance Flux 

Source 

(1000 km 2 ) 

(Gt a 1 ) 


159 

63.4 

From DEM by Liu et al. [1999] and surface 
balance compilation of Vaughan et al. [1999] 

175 

69 

From DEM by Bamber and Bindschadler [1997] 
and surface balance compilation of Vaughan et al. 
[1999] 

159 ± 1 

63.9 ±6 

[Rignot y 1998] 


76 

[Bentley and GiovinettOy 1991] 

(Arithmetic mean of estimates from Crabtree and 
Doake [1982] and Lindstrom and Hughes [1984]) 

182 

65. 9± 5 

[Lindstrom and Hughes t 1984] 

214 ±20 

86± 30 

[Crabtree and Doake y 1982] 


the southern lobe , feeding the first through a neck less than 
100 km across. A delineation of the catchment basins of 
the 70 largest glaciers in Antarctica, by a similar method, 
showed this configuration is unusual [Vaughan and 
Bamber , 1998]. Generally, ice-drainage basins which drain 
through glaciers or ice streams, are uniformly convergent. 
Only the Ice Stream C basin has a similar “necked” shape 
[. Joughinetal , 1999]. It is possible that the existence ofthe 
southern lobe of the Pine Island Glacier basin indicates 
unsteady conditions in the basin, with this lobe currently 
being transferred between catchment basins. Alternatively, 
it may simply reflect unusual bed morphology. 

2.3 Mass input and balance flux 

Overlaying the basins for Pine Island Glacier derived in 
Section 2.2, on a grid representing the mean surface 
balance over Antarctica [Vaughan et al y 1999], we have 
estimated the total rate of snow accumulation in the Pine 
Island Glacier basin (Table l). This is the amount of ice 
that must leave the basin for mass balance to be maintained 
and is termed the balance flux. The aggregate of the three 
most recent estimates gives a balance flux for Pine Island 
Glacier of (66 ± 4) Gt a 1 . The uncertainty is derived from 
the spread of the results, but is consistent with the 
uncertainties in area (± 4%) and accumulation (± 5%) 
[ Vaughan et al. , 1 999] . 

A similar process has been used to find the balance 
fluxes of the other major glaciers of Antarctica, and while 
many glaciers are fed by basins with larger areas, only one 


balance flux exceeds that of Pine Island Glacier. Totten 
Glacier, East Antarctica has a balance flux of around 75 Gt 
a' 1 [Vaughan and Bamber , 1998]. Outside Antarctica, the 
most active glacier is Jacobshavns Isbrae, Greenland which 
supports about half this flux [Bindschadler, 1984]. 

2,4 Glacier tributaries 

Two techniques employing satellite data and yielding 
wide coverage allow us to identify areas where ice-flow is 
concentrated within the interior of the ice-drainage basin. 
The pattern that emerges is one of great complexity, with 
many tributaries coalescing to form the main glacier. 

2 . 4.1 Satellite Altimetry. We can derive some 
understanding of the distribution of ice flow in the interior 
of the basin using the OEMs discussed in Section 2. 1 . The 
method first calculates the flow-direction for each cell. It 
then assigns to each cell a numerical value corresponding 
to the number of other cells whose accumulation will 
eventually flow through it. The technique is known as 
flow-accumulation and a grey-scale representation of this 
flow-accumulation grid (Figure 5) gives an indication of 
where flow is more convergent within the basin 

Figure 5 shows a system of tributaries which merge 
about 100 km above the grounding line to form the single 
unit of flow which is Pine Island Glacier. These tributaries 
are identifiable several hundreds of kilometers inland, 
much further inland than the point where the ice enters a 
more confined channel, previously suggested to be the start 
of channelized ice flow [Lucchitta et al. , 1 995; Lucchitta et 
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Fig. 5. Map of area around the basin drained by Pine Island Glacier showing flow-accumulation 
derived from 5-km resolution surface elevation grid. Grid cells are shaded such that cells fed by 
many others are darker than those fed by only a few. The darker areas thus represent areas into 
which the flow is channelled The numbered features are the tributaries as identified by Stenoien 
f 1 998] Stenoien’s tributaries l and 10 arc not resolved on this representation. 



al, 1 994]. This set of tributaries is almost the same as that 
determined by Stenoien [1998] from interferometric 
Synthetic Aperture Radar (SAR) images (see figure 6.14 
of Stenoien [1998]) and in Figure 5 they are numbered 
using Stenoien’s designation. One tributary 7 (5) drains 
much of the southern lobe of the catchment basin 
described in Section 2.2 Its presence is perhaps not 
surprising given the narrow neck where the southern lobe 
joins the rest of the drainage basin. 

Stenoien [1998] suggested that similar patterns of 
tributaries have not been seen elsewhere, and that they are 
perhaps unique to Pine Island Glacier. There is, however, 
evidence elsewhere in Antarctica for tributary systems in 
other basins. Bright margins in ERS-1 SAR data have 
shown that Evans Ice Stream forms at the confluence of at 
least five tributaries [Jonas and Vaughan , 1 996]. Radarsat 


data shows that Recovery Glacier has two tributaries w hich 
extend hundreds of kilometers inland [ Jezek , et al.. 1998]. 
Flowlines in Landsat imagery- show that Institute Ice 
Stream also has several tributaries [Mantripp et al. , 1 996]. 
Finally, Joughin et al [1999] show a tributary system 
feeding ice streams on the Siple Coast. 

The presence of several tributaries in the Interior 
drainage system may imply that Pine Island Glacier is 
unlikely to respond dramatically to changes in one locality. 
For example, if “water-piracy” [Alley et al, 1994] or a 
reduced supply of basal till were to shut-off one of the 
tributaries, the others would probably be unaffected and 
the flux in Pine Island Glacier may suffer little change. 

2*4.2 Interferometric SAR. Goldstein et al [1993] 
used SAR data from the ERS-1 satellite to construct 
interferometric SAR (InSAR) images of ice flow in 
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Antarctica. These showed ice movement only along the 
line of sight to the satellite, but the method has since been 
refined to produce a 2-dimensional velocity-field for parts 
of the interior of the Pine Island Glacier basin 
[Stenoien, 1998]. The procedure yields an understanding 
that is more quantitative than that from the satellite 
altimetry presented above, although coverage of the Pine 
Island Glacier basin is less complete. Stenoien ’s [1998] 2- 
dimensional velocity field covers much of the northern 
part of the drainage basin (see figures 5.14 and 5.15 of 
Stenoien [1998]). The data cover die upstream regions of 
many, of the tributaries (2, 4 and 6 in Figure 5) and the 
middle regions of two which originate in the southern lobe 
of the basin (3 and 5 in Figure 5). There are no points of 
ground control in this i-ea and Stenoien derived an 
absolute velocity field by assuming that a saddle on the ice 
divide had zero velocity. Ice speed away from the 
tributaries is low (0-50 m a 1 ) but increases within them to 
more than 150 m a' 1 upstream of the confluence of 
tributaries 2, 3 and 5 (Figure 5). 

Of interest is Stenoien’ s observation that none of the 
tributaries for which data are available shows a rapid 
increase in ice speed, but rather a gradual increase down 
the length of each tributary. A similar pattern has also been 
observed on other West Antarctic ice streams (Joughin et 
al , 1 999) The lack of a sudden velocity increase, suggests 
notions of a bi-stable state of glacier-flow i.e. fast or slow, 
may be unrealistic, but rather that a progressive response 
to changing boundary conditions is possible. 

Taken together, InSAR and flow-accumulation show 
that the interior of the Pine Island Glacier basin is complex 
with around 10 tributary ice streams coalescing to form a 
single glacier. None of these tributaries appears to have a 
well-defined region of rapid velocity increase and what 
controls their location and longevity remains to be 
determined. Radar data presented in Section 2.6 suggest 
that the control may be through basal conditions. However, 
at present even sophisticated thermomechanical models of 
the area fail to reproduce this complex flow pattern [e.g. 
Payne , 1999]. 

2,5 Subglacial topography 

The bed topography of the West Antarctic ice sheet was 
first mapped in detail using a combination of traverse data, 
airborne sounding data and TWERLE balloon altimetry 
[Jankowski and Drewry, 1981], While this study clearly 
delineated the major subglacial features of the area, the 
availability of new data prompts us to repeat the exercise. 

Figure 6 shows a new compilation of bed topography 
beneath the grounded portion of the Pine Island Glacier 
basin. To create a grid of ice thickness we used (Figure 6a) 


traverse data [Bentley and Ostenso, 1961; Behrendt, 1964; 
Bentley and Change 1971], airborne radar data [Jankowski 
and Drewry, 1981; British Antarctic Survey unpublished 
data] and rock outcrops, which were used as an isopleth of 
zero ice thickness. The resulting grid of ice thickness was 
subtracted from the ERS-1 derived DEM of surface 
elevation described in Section 2.1 to produce bed 
topography (Figure 6b). 

The bed topography (Figure 6b) shows clearly the main 
features identified in earlier compilations: the Bentley 
SubglaciafTrench; the Byrd Subglacial Basin, which here 
reaches almost 2000 m below sea level; and between these 
depressions” the “sinuous ridge” described by Jankowski 
and Drewry [1981]. The Ellsworth Subglacial Highlands, 
are also well-defined. Despite significantly improved data 
coverage ; n the Pine Island Glacier basin. Figure 6b shows 
no new substantial features except a trough 1 000 m below 
sea level, in which Pine Island Glacier and its main 
tributaries flow. 


2,6 Driving stress 

The driving stress in an ice sheet is calculated from the 
surface slope and ice thickness according to a simple 
relation [Paterson, 1994; page 241]. Here we have 
calculated the driving stress for the region (Figure 7) using 
the ERS-1 -derived DEM and the ice thickness grid 
described above. The calculated driving stress is negligible 
near the ice divide where the surface slopes are low; 
intermediate (50-110 kPa) on the slow-moving areas 
between the ice divides and the tributary glaciers; and low 
(<50 kPa) on the tributary glaciers, but rises to >1 10 kPa 
along the main trunk of Pine Island Glacier. 

An airborne sortie was flown from the inactive Siple 
Station (75° 54’ S 84° 30' W) to the ice front of Pine Island 
Glacier in 1998 (the flight-track is shown in Figure 6a). It 
covered much of the main tributary of Pine Island Glacier 
(that formed by numbers 2, 4, 6 and 8 in Figure 5). Ice- 
penetrating radar data from this sortie show that the margin 
of this main tributary is marked by a downward step in the 
bed elevation and a change to a smoother ice-base 
reflection which has an “ice shelf-like” character (Figure 
8). This change in reflection character is believed to 
indicate a transition from a frozen bed (rough) to one 
which is at the pressure-melting point (smooth). Although 
the flight track does not follow an ice flow line exactly, the 
driving stresses derived from the along-track data compare 
well with those derived from the gridded datasets shown in 
Figure 7. The driving stress calculated from the along- 
track topography has four distinct zones (Figure 8); the 
interior of the basin (50-75 kPa), the main tributary glacier 
(around 30 kPa), the main trunk of Pine Island Glacier 
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Fig. 6a. Map of the area around the ice-drainage basin of Pine Island Glacier -showing 
measurements of ice thickness from airborne survey and oversnow traverses. Thin g y 
indicate unsuccessful sounding by airborne survey. 


(over 100 kPa) and the floating ice (less than 10 kPa). The 
marked change in bed roughness across the tributary 
margin may also indicate that the location of the tributary 
and underlying geological constraints are closely related. 
The low driving stress suggests that the tributary flows 
over a well-lubricated bed. 

The pattern of driving stresses suggests that the Pine 
Island Glacier basin is dynamically different from the 
idealised ice-stream basin. Much of the basin comprises a 
slow-moving ice sheet which may be cold-based, as 
suggested by the character of the radar reflection (Figure 
8). This ice sheet feeds a number of wet-based (Figure 8), 
lubricated tributaries with relatively low driving stresses 
(around 30 kPa). These merge to form Pine Island Glacier, 
which has a much higher driving stress (>100 kPa), more 
akin to an East Antarctic outlet glacier, than some of the 
West Antarctic ice streams [ Bentley , 1987], 


2.7 Surface elevation change 

ERS- 1 satellite altimetry data for the period 1 992- 1 996 
were analysed for evidence of surface-elevation change 
[Wingham et a/,1998] The data covered most of the 
interior of the Antarctic Ice Sheet north of 82°S. The 
analysis showed only one region of spatially-coherent 
surface-elevation change. Thinning at a mean rate of 11 .7 
± 1 .0 cm per year was indicated in the Pine Island Glacier- 
Thwaites Glacier basin. Wingham et a! [1998] indicated 
that the change was centered and most significant over the 
Thwaites Glacier basin (see their Figure 2), rather than the 
Pine Island Glacier basin, but the trend did appear to 
extend across both. The simplest interpretation is that the 
surface lowering resulted from a change in surface mass 
balance. Alternatively, a change in the glacier flux due to 
increased discharge or grounding-line retreat might also be 
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Fig. 6b. Contour map of bed elevation (meters above sea level). A grid of ice thickness was 
calculated using the measurements shown in Figure 6a. This grid was subtracted from a surface- 
elevation grid [Bamber and Bindschadler , 1997] to produce a grid of bed elevation. 



the cause. However, in either case, the shortness of the 
observation period gives little indication of future 
behavior. It is hoped that more detailed analysis of the 
ERS-1 altimetry will refine the pattern of change. NASA's 
Geoscience Laser Altimeter System (GL AS), scheduled for 
launch in 2001, will allow similar measurements to be 
made even in the coastal margin of Antarctica. 

3. THE GLACIER 

Until the discovery of the network of ice tributaries in 
the basin [Stenoien, 1 998; Section 2.4], Pine Island Glacier 
was generally considered to extend only around 70 km 
above the grounding line to where the ice is first 
channelled into parallel flow (Figure 9). This may still be 
still a useful definition, since it draws some distinction 
between the tributaries and the main trunk of the glacier 


and approximately marks the increase in driving stress 
mentioned in Section 2.6. Thus defined, Pine Island 
Glacier is bounded to the north by nunataks in the Hudson 
Mountains and to the south by slow-moving ice sheet. 

3.1 Surface features 

Surface features on Pine Island Glacier revealed by 
Landsat and SAR imagery have been shown by various 
authors and are reproduced in Figure 9. Flowlines of the 
type discussed by Whillam and Merry / [1993] show 
considerable convergence at the head of the glacier, around 
the zone of arcuate “crevasses” revealed by ERS-SAR 
images (shown in Figure 9 and in greater detail by 
Lucchitta et al. [ 1 995]). These presumably mark a zone of 
longitudinal extension. Lucchitta et al [1995] noted that 
these “crevasses” had not been previously described and 
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Fig. 7. Map of driving stress for the grounded ice sheet in Pine Island Glacier and neighbouring 
ice-drainage basins derived from surface elevation (Figure 3) and bed elevation (Figure 6b). Light- 
grey shading denotes driving stresses in the range 0-50 kPa, mid-grey 50-1 10 kPa, and dark-gre> 
greater than 1 10 kPa. 


are not shown by visible images, and concluded that they 
may be covered by a layer of snow. 

Below the zone of convergence and arcuate crevasses, 
flowlines on the glacier are roughly parallel. The surface 
has smooth, long-wavelength (a few kilometers) 
undulations that are typical of fast-moving glaciers. 

3.2 Velocity 

Lucchitta and others [Lucchitta et a /., 1995; Lucchitta 
et al. , 1994; Ferrigno et al. , 1993] have measured the ice 
velocity on Pine Island Glacier using sequential SAR 
images acquired by the ERS- 1 satellite. They found that on 
the main trunk of the grounded portion, the center-line 
speed ranged from about 1 km a' 1 near the arcuate 
crevasses, to 1.5 km a' 1 at the grounding line identified by 
Crabtree and Doake [1982], The flow speed then rose 


rapidly to 2.5 km a' 1 between that grounding line and the 
one identified by Rignot [1998], and then remained 
approximately constant to the ice front. Their velocity 
measurements were generally higher than earlier estimates 
[Kellogg and Kellogg, 1987; Lindstrom and Tyler, 1984; 
Crabtree and Doake, 1982. Williams et al., 1982] although 
the data are not adequate to determine if there has been any 
acceleration. 

3.3 Grounding line 

The grounding line of Pine Island Glacier is not easily 
discemable on either the Landsat or ERS-1 SAR imagery 
in Figure 9, but its position was determined by hydrostatic 
calculations based on airborne data [Crabtree and Doake, 
1982] (Figure 9). This hydrostatic condition downstream 
of this grounding line was, however, not entirely clear. 
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Fig. 8. Driving stress, topography and radar data collected during airborne survey flight-track 
identified in Figure 6a. a) Driving stress calculated using 10-km smoothing along the flight-track. 
The driving stress is calculated to be positive in the direction of Pine Island Glacier ice front and 
is negative where the ice flows into a different basin, b) Ice-surface and bottom topography along 
the flight-track with ice-flow features marked, c) Section of radar data, showing the difference in 
character between the ice-bottom return from the interior drainage basin (rough) and fiom the 
tributary glacier (smooth). This difference probably reflects the difference between a frozen bed 
and one that is at the pressure-melting point. 


Thomas [1984] argued that a zone of partial grounding 
might exist for 30 km downstream. This downstream 
position has now been confirmed, using InSAR to detect 
the limit of tidal flexing [Rignot, 1998]. That analysis also 
showed that the limit of flexure extends seaward (around 
1 5 km) near the centre of the glacier with a re-entrant on 
either side (Figure 9). This pattern is possibly an indication 


of a bedrock obstruction similar to that known to underlie 
Rutford Ice Stream [, Stephenson , 1984; Doake et al, this 
volume] or it might be solely due to the glacier being 
thicker close to its center line. 

The position of the limit of flexure was measured at 
five epochs between 1992 and 1996, which showed that it 
moved during this period [ Rignot , 1998]. The simplest 
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Fig, 9a. Landsat 1 sub-scene of Pine Island Glacier acquired on January 24 } 1 973 (path 246, row 
1 1 4). marks the position of the raft of ice discussed in Section 4. 1 . “U" indicates plumes of 
^ ice thickness undulations formed close to the grounding line and dissipating towards the ice front 

(discussed in Section 4.1). The dotted line marks the grounding line identified by Crabtree and 
Doake [1982] and the black line indicates the limit of tidal flexing determined by Rignot [1998] 
for 21 January' 1996. 

Fig. 9b. Mosaic of two ERS-I SAR images (orbit 3174, frames 5193 and 521 1) acquired 
December 4, 1992 showing the same area. “A'’ marks the set of arcuate crevasses identified by 
Lucchitta et al [1995]. 
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interpretation is that between 1992 and 1994 there was a 
retreat in the position of the grounding line in the center of 
the glacier that averaged 1 .2 ± 0.3 km a* 1 . The pattern of 
change within the re-entrant parts is not so clear, though 
Rignot calculated that it could be caused by a thinning rate 
of 3.5 ± 0.9 m a' 1 . This is a small fraction of the basal melt 
rates in this area (see Section 4.2) and thus could result 
from a relatively small change in the oceanographic 
conditions. Alternatively, it could also be caused by a 
thinning of the glacier upstream of the grounding line. 

4. THE ICE SHELF 

Pine Island Glacier debouches into an ice shelf 
comprising the ~80 km floating portion of the glacier plus 
the slow-moving floating ice sheet that surrounds it. The 
floating portion of Pine Island Glacier is easily identified 
by plentiful flowlines generated near the grounding line, 
which continue to the ice front showing little divergence. 

4.1 Surface features 

Kellogg etai [1985] found that dense (650 kg 
nr 3 ) “well-sintered” fim predominated at the surface of 
the floating portion of Pine Island Glacier, close to the ice 
front. They interpreted this as resulting from strong 
katabatic winds, producing net sublimation from the ice 
surface. The extent, persistence or magnitude of this 
negative net surface mass balance is, however, not well- 
established. 

Landsat and SAR images of the floating portion of Pine 
Island Glacier (Figure 9) show crevasses formed at the 
grounding line, moving in plumes to the ice front. Much of 
the southern side of the floating portion of the glacier is 
covered by periodic transverse surface undulations visible 
on the Landsat images. These also form in plumes 
emanating from the grounding line and dissipating towards 
the ice front. Airborne radar-sounding data show that these 
undulations have an amplitude of around 20 m, a 
horizontal wavelength of around 2.5 km and are 
hydrostatically compensated by ice thickness changes of 
190 m. The ice flow velocity measured on this section of 
ice shelf was 2.3 - 2.6 km a 1 [Lucchitta et al. t 1997] which 
suggests that one undulation is produced each year, 
although the mechanism that causes them is uncertain. 

A former feature of the floating portion of Pine Island 
Glacier, not previously discussed but clearly visible in 
Figure 9, was a raft of thicker ice embedded in the ice 
shelf. In later images, the raft was seen to have been 
advected downstream at approximately the same speed as 
the ice-shelf flow. The origin and significance of this raft 


is unclear and it is no longer available for investigation as 
it would have calved from the ice shelf in the late- 1980s. 
Similar features have been noted in grounded ice streams 
[e.g., Whillans et al, 1993] and possibly in floating ice 
shelves [Cassasa etal, 1991]. 

4.2 Basal melting 

In early 1994, the research ship Nathaniel B. Palmer 
entered Pine Island Bay and conducted an oceanographic 
survey of the area which has led to three studies revealing 
a most unusual oceanographic regime: 

• Jacobs et al [ 1 996] used oceanographic measurements 
and a “salt-box” calculation to show that the mean basal 
melt rate beneath the floating portion of Pine Island 
Glacier was 10-12 m a 1 , five times the highest rate 
previously published (on George VI Ice Shelf [Bishop 
and Walton , 1981]). They concluded that these high 
melt-rates were driven by relatively warm Circumpolar 
Deep Water (CDW) flooding this portion of the 
continental shelf, combined with the deep draft of Pine 
Island Glacier. 

• Jenkins et al. [ 1 997] determ ined that the position of the 
ice-shelf front showed no persistent trend in the period 
1973-1994 (in Figure 10 we extend this series to 1966- 
1998). Combining this with a flux calculation across the 
grounding line and with the assumption that it is a 
steady-state system, they calculated a mean basal melt 
rate over the ice shelf of 12 ± 3 m a* 1 . 

• Hellmer et al [ 1998] used analyses of dissolved oxygen 
and oxygen isotopes to confirm the strong melt-water 
signal in the outflow and applied a thermohaline model. 
This suggested that in some areas, the basal melt-rate is 
twice the mean value and that temporal variations in the 
temperature of the inflowing CDW could cause 
substantial changes in the basal melt-rate. 

Rapid melting from beneath the floating portion of Pine 
Island Glacier jwas_ later confirmed using satellite 
measurements of mass balance [Rignot, 1998] giving a 
mean melt-rate of (24±4) m a* 1 increasing to (50± 1 0) m a ‘ 
close to the grounding line. These melt-rates are larger 
than those calculated by Jenkins et al [1997] due to a new 
position for the grounding line constrained by InSAR 
observations of tidal flexing. Using this grounding line 
position, we increase the estimate of mean melt-rate 
produced by Jenkins et al. to around 17 m a 1 , and that 
from Jacobs et al. by a similar amount. 

These observations of a thick ice shelf coming in 
contact with relatively warm Circumpolar Deep Water has 
led to speculation that ice-ocean interactions in Pine Island 
Bay may be similar to those prevalent during the Last 
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Fig, 10. Map of selected ice-front positions for Pine Island Glacier, between 1966 and 1998 
overlaid on an excerpt of a sketch-map drawn from aerial photography collected in 1 966 (USGS, 
1 993 ; original map prepared in 1 967). Sources: 1 966, Aerial photography (USGS, 1 993); January' 
24,1973, Landsat 1 (path 246, row 1 14); January 15, 1982, Landsat image; February 9, 1992, 
ERS-l SAR image; December 4, 1992, ERS-I SAR images (orbit 3 1 74, frames 5193 and 5211); 
15 March 1994, ERS-l SAR image; February, 1996, ERS-I SAR image; February 13, 1998, 
flight-track of BAS airborne survey aircraft which flew along the ice front. 
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Glacial Maximum when much of the West Antarctic ice 
sheet probably extended out to the edge of the continental 
shelf and came into contact with similarly warm water [e.g. 
Jenkins et al , 1 997]. In Pine Island Bay we now know that 
such conditions generate high sub-ice shelf melt-rates and 
we infer that these rates may have also been prevalent 
during glacial periods, severely limiting the size of any ice 
shelves. In addition, the results highlight the importance of 
oceanographic conditions as a significant control on the 
present and future configuration of the West Antarctic ice 
sheet. 

4.3 Ice front stability 

The ice-front position of Pine Island Glacier has been 
reconstructed by several authors using several sources of 
data and covering various periods. Figure 10 shows the 
position of the ice front sporadically since 1966. The 


pattern suggests a retreat of around 10 km between 1966 
and 1973, followed by a period of general stability, and 
readvance in recent years. However, this simple 
interpretation should be qualified on two counts, a) The 
presence of an iceberg just off the ice front on the sketch 
map drawn from 1966 aerial photography. This indicates 
that a significant calving had recently occurred, and that 
prior to this the ice front was even further advanced than 
the most extreme position shown. The same is true for the 
1973 image, b) The ice velocity at the front (>2.5 km a*’) 
allows for fluctuations within the intervals between the 
observations, large enough to exceed the extreme positions 
shown in Figure 10. However, it is interesting that the ice 
front position appears to be relatively stable, despite the 
high ice velocity. Further monitoring would be required to 
confirm whether or not this is the case. 

Kellogg and Kellogg [1987] have suggested that the 
Thwaites Iceberg Tongue was not actually formed from 
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calving of Thwaites Glacier Tongue but might have calved 
from Pine Island Glacier. This is, however, unlikely, as a 
prolusion of surface transverse lineations on the iceberg 
tongue matched well with the transverse lineations on 
Thwaites Glacier, but did not match the longitudinal 
lineations that predominate on Pine Island Glacier 
[Ferrigno et al, 1993]. 

In summary, there are insufficient data available to 
discern a decadal trend in the ice-front position of Pine 
Island Glacier, although it is probable that some retreat 
(-10 km) has occurred in the last 30 years. (Longer-term 
retreat of the ice shelf front is discussed in Section 5). 

5. THE MARINE ENVIRONMENT 
5. 1 Retreat of ice in Pine Island Bay 

Seabed sediments provide information on glacier retreat 
in Pine Island Bay. Results from four cores within Pine 
Island Bay and from 1 9 cores on the outer continental shelf 
and eastern Amundsen Sea have been presented [Anderson 
and Myers, 1 98 1 ; Kellogg and Kellogg, 1987; Kellogg et 
al, 1987], but their value is limited because they contain 
very little material suitable for radiocarbon dating. The 
radiocarbon dates that do exist, only poorly constrain the 
retreat of the ice sheet in Pine Island Bay to the last few 
millennia. Prior to this, the ice sheet may have occupied 
the entire Pine Island Bay, and perhaps butted against the 
Thwaites Glacier Tongue to form an extensive ice shelf or 
ice sheet. 

Kellogg and Kellogg [1987] suggested the retreat was 
very recent, with Pine Island Bay being filled with 
grounded ice only 100 years ago, but this interpretation 
relied heavily on an ongoing retreat rate of 0.8 km per 
year inferred from aerial photography acquired in 1966 
and Landsat imagery acquired in 1973. The variable 
position of the ice front shown in Figure 10 now casts 
doubt on such an extrapolation. 

5.2 Ice Extent at the Last Glacial Maximum 

The sediment cores collected on the outer continental 
shelf, near 1 10'W [Anderson and Myers, 1981], and north 
of Thurston Island between 100 - W and 102‘W [Kellogg 
and Kellogg 1 987; Kellogg et al . , 1 987] show a thin (0- 1 5 
cm) upper layer of sandy mud, probably of Holocene age, 
containing common planktonic and calcareous benthic 
foraminifera. Diatoms are relatively rare in this layer 
despite high abundances in the surface water. A compact, 
poorly-sorted diamicton underlies the sandy mud. This 
diamicton is generally more than 2.3 m thick, and probably 


represents deposition beneath grounded ice that extended 
to the continental shelf break. 

If we assume that the diamicton is a remnant of 
subglacial basal till, then grounded ice probably remained 
over most of the Amundsen Sea continental shelf until 
relatively late in the Holocene. The postglacial sediments 
on the outer shelf are much thinner in Pine Island Bay than 
on the outer shelf in the Ross Sea [Kellogg et al, 1979; 
Domacket al, 1999; Licht et al, 1999; Shipp et al ., 1999] 
which, given similar deposition rates, might suggest an 
earlier deglaciation in the Ross Sea than in the Amundsen 
Sea. Finally, if we assume that post-glacial deposition rates 
were similar to current measured rates (e.g., 10-35 cm a 1 
beyond coastal Alaskan glaciers [Molnia and Carlson, 
1999] and ~10 cm a' 1 in Antarctic Peninsula fjords 
[Domack and McClennen, 1996]), then we can conclude 
that deglaciation occurred, at most, a few thousand years 
ago. 

6. MASS BALANCE 

Comparisons of balance flux (Table 1) and grounding 
line flux (Table 2) for Pine Island Glacier show a 
progression toward reduced uncertainty, but significant 
uncertainty still remains. Our preferred estimate of overall 
mass balance is -2.4 ± 4 Gt a 1 (the difference between the 
balance flux determined in this study and the grounding 
line flux calculated by Rignot [1998]). This value is not 
significantly different from zero. Given a catchment basin 
area of -170 000 km 2 , this would be equivalent to a 
lowering of surface elevation in the range 1.5 - 3 cm a 1 , 
depending on the density of the layers being lost. 

Rignot [1998] determined the ice thickness by inverting 
the ice-surface elevation at the limit of flexing using a 
hydrostatic condition but there is evidence that the limit of 
flexing is often many kilometers upstream of the 
hydrostatic limit. One example is on Rutford Ice Stream 
where along the center-line the limit of flexing is 2 km 
upstream of the hydrostatic point and the surface is 50 m 
above the hydrostatic condition [ Vaughan , 1994; Smith, 
1991]. If a similar situation applies on Pine Island Glacier 
then the ice flux across the grounding line may have been 
over-estimated by Rignot. However, as at present we 
cannot assess the likelihood of this possible error, we use 
Rignot’s grounding line flux for our preferred mass 
balance estimate. 

The mass balance within two regions of the drainage 
basin has also been calculated (see figure 6.8 and table 6. 1 
of Stenoien, [1998]). The north-eastern part of the basin 
has a positive mass balance (6.4±3.7 Gt a '), whereas in the 
region just north of the narrow neck it is negative 
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Table 2. Estimates of grounding-line flux for Pine Island Glacier. 


, Mass (Gt a* 1 ) 

Source 

68.4 ±2 

[Rignot, 1998] 

>56 ±6 

[Jenkins et al , 1997] 

70 

[Lucchitta et al , 1995] 

25.5 ±5 

[Lindstrom and Hughes , 1984] 
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(-7/7±4.7 Gt a* 1 ). These values were calculated using a 
hand-drawn map of mean surface mass balance, but 
repeating the calculation using an updated map of surface 
mass balance [Vaughan et a/., 1999], gave the same 
answer. Distributed evenly across the areas, these 
imbalances represent changes in surface elevation of 
(23±14) cm a 1 and (-51±3l) cm a' 1 , respectively. 

There is thus a contradiction between mass balance 
calculations and measurements of change in surface 
elevation [Wingham et a/., 1998; see Section 2.7] which 
seems to imply one of three possibilities; a) substantial 
changes in the density-depth relation in the snow in the 
period 1992-1996, b) unusually low precipitation 
accumulation in the period 1992-1996 or, c) one, or both 
of the analyses is substantially in error. 

7. DISCUSSION AND CONCLUSIONS 

The Pine Island Glacier ice-drainage basin is an area of 
particular interest and in some ways, may be unique. The 
shape of the basin is unusual. It comprises two lobes joined 
by a narrow neck less than 100 m across. The southern 
lobe has a higher ice-surface elevation and does not appear 
to correlate with any significant bed feature. Surface 
elevation in the ice-drainage basin dropped over four years 
[Wingham et ai, 1998], though the reasons for this are 
uncertain. Precipitation rates are high and the basin has the 
second-highest balance flux of any extant ice stream or 
glacier (66±4 Gt a' 1 ). Although there are indications that 
the ice sheet may be out of balance locally [Stenoien, 
1998], mass balance calculations of the basin as a whole 
show that there is currently no measurable, significant 
imbalance. The tributaries which drain the basin flow at 
intermediate speeds (50- 1 50 m a' 1 ) and show no regions of 
marked velocity increase. They coalesce to form Pine 
Island Glacier which has high driving stresses (> 1 00 kPa) 
similar to East Antarctic outlet glaciers, but also shear 
margins bordering slow-moving ice, often typical of West 
Antarctic ice streams. The bed of the slow-moving ice in 


the drainage basin appears to be frozen. In the tributaries, 
the bed is at the pressure-melting point and is well- 
lubricated, presumably associated with the reduction 
observed in the driving stress. This thawed bed continues 
into Pine Island Glacier. The glacier flows over a complex 
grounding zone where Rignot [1998] measured a 
grounding line retreat along the center of the glacier of a 
few kilometers over a period of a few years. Once afloat, 
the base of the glacier comes into contact with warm 
Cicumpolar Deep Water which generates basal melt-rates 
of 10-20 m a' 1 or more [Jacobs et ah , 1996; Jenkins et al , 

1 997; Hellmer et al. , 1998 \Rignot. 1998], the highest rates 
yet measured beneath any ice shelf The front of the ice 
sheet or shelf in Pine Island Bay may have retreated 
significantly in the last few millennia [Kellog and Kellogg 
1987], but over the last few decades it appears to have 
been more stable, shifting back and forth no more than a 
few kilometers. 

These observations summarize what is known about 
Pine Island Glacier and its drainage basin. Several of them 
indicate evidence of change in this part of the West 
Antarctic ice sheet on time scales varying from millennia 
to a few years. However, at present we are cautious as to 
the long-term significance of these changes. Although we 
note that the unusual basin shape could indicate ongoing 
transfer of the southern lobe between catchments, there is 
no evidence to support this and the shape may simply be 
reflecting the bed topography in some way. We are 
uncertain as to the whether the observed surface-elevation 
change is the result of changing precipitation or else 
changing ice flow. In addition, the apparent local 
imbalances do not seem to concur with the changes in 
surface elevation. We do not understand yet if the recent 
grounding line retreat has been associated with a 
significant change in the force-balance of the glacier. 
Although the ice-ocean interaction in Pine Island Bay 
appears to be unusually dynamic, we cannot be sure of the 
past or future nature of the intrusion of warm water which 
is responsible. Finally, more radio-carbon dating and better 
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discrimination between ice-rafted, sub-ice shelf and sub- 
ice sheet sediments are required to improve our confidence 
in the indications of ice-front retreat over the past 
millennia. Hence, as yet, none of the observed changes 
make a strong case for ongoing basin-scale ice-sheet 
change or readjustment and certainly, none suggest that the 
ice sheet in this area has entered a phase of significant 
collapse or retreat. 

At a time when the paradigm of marine-ice-sheet 
instability is being questioned, our observations of change, 
perhaps appear ambiguous and inconclusive. Interpreting 
them as precursors of collapse would clearly be 
unjustified. They are certainly not yet sufficient to 
rigorously test the various theoretical models and ideas of 
ice sheet stability and collapse which have been discussed 
over the years [e.g. Weertman , 1974; Mercer , 1978; 
Hughes , 1980; Fas took, 1984; MacAyeal , 1992; 

Hindmarsh , 1993; Bentley , 1998]. This is the case, even 
though Pine Island Glacier has occupied a central position 
in discussions regarding the stability of the West Antarctic 
ice sheet. It can be argued that research has made relatively 
poor progress in the Pine Island Glacier region, 
particularly when compared to the Siple Coast ice streams 
or those which drain into the Ronne Ice Shelf. This 
shortcoming suggests that we should prepare field and 
remote sensing experiments that will allow us to determine 
the causes of change in the ice sheet. Only such an 
understanding will lead us to a sound foundation for 
predicting future behaviour. 

Acknowledgments . We thank Howard Conway and Robert 
Bindschadler for constructive reviews. SSJ acknowledges support 
of the NASA Polar Research Program for the LDEO contribution 
(#6034). 

REFERENCES 

Alberts, F.G., Geographic names of the Antarctic. National 
Science Foundation, Washington. 1981 
Alley, R.B., S. Anandakrishnan, C.R. Bentley and N. Lord, A 
water-piracy hypothesis for the stagnation of Ice Stream C, 
Antarctica, A nn. Glaciol , 20, 187-194, 1994 
Anderson, J.B. and M.C. Myers, USGSC Glacier Deep Freeze 81 
Expedition to the Amundsen Sea and Bransfield Strait, 
Antarct. J U.S. , 16, 5, 1981 

Armstrong, T., B. Roberts and C W.M. Sw ithinbank, Illustrated 
glossary' of snow and ice. Scott Polar Research Institute, 
Cambridge, 1973. 

Bamber, J.L. and R.A. Bindschadler, An improved elevation 
dataset for climate and ice-sheet modelling: validation with 
satellite imagery, Ann. Glaciol, 25, 439-444, 1997. 
Behrendt, J.C., Distribution of narrow-width magnetic anomalies 
in Antarctica, Science , 144, 3621, 993-994, 1964. 


Bentley, C.R., Antarctic ice streams: a review, J Geophys. Res., 
92, 8843-8858, 1987. 

Bentley, C.R., Rapid sea-level rise from a West Antarctic ice- 
sheet collapse: a short-term perspective, J. Glaciol., 44, 1 57- 
163, 1998. 

Bentley, C.R. and N.A. Ostenso, Glacial and subglacial 
topography of West Antarctica. J. Glaciol , 3, 882-91 1 . 1961. 

Bentley, C.R. and F.K. Chang, Geophysical exploration in Marie 
Byrd Land, Antarctica, in Antarctic Snow and Ice Studies 2, 
Antarct. Res. Ser. vol 16, edited by A.P Crary, pp 1-38, 
AGU, Washington, D.C., 1971. 

Bindschadler, R., Jakobshavn Glacier drainage basin: a balance 
assessment, J. Geophys. Res., 89, 2066-2072, 1984. 

Bishop, J.F. and J.L.W. Walton, Bottom melting under George 
VI Ice Sheet, Antarctica, J. Glaciol , 27, 429-447, 1981. 

Bromwich, D. Snowfall in high southern latitudes. Rev 
Geophys., 26, 149-168, 1988. 

Cassasa, G., K.C. Jezek,, J. Turner and LM. Whillans, Relict flow 
stripes on the Ross Ice Shelf. Ann. Glaciol , 15, 132-138, 
199! 

ConnoIIey, W.M. and J.C. King, A modelling and observational 
study of East Antarctic surface mass balance, J. Geophys 
Res., 101, 1335-1343, 1996. 

Crabtree, R.D. and C.S.M. Doake, Pine Island Glacier and its 
drainage basin: Results from radio-echo sounding, Ann. 
Glaciol, 3, 65-70, 1982. 

Doake, C.S.M., H.F.J. Coit, A, Jenkins, K. Makinson, K.W 
Nicholls, C. Nath, A.M. Smith and D.G. Vaughan, Rutford 
Ice Stream, Antarctica, this volume . 

Domack, E. W. and C. E. McClennen, Accumulation of glacial 
marine sediments in fjords of the Antarctic Peninsula and 
their use as late Holocene palaeoenvironmental indicators, in 
Foundations for Ecological Research West of the Antarctic 
Peninsula, Antarct. Res. Ser. vol 70, edited by R.M. Ross et 
ai., pp 135-154, AGU, Washington, D C., 1996. 

Domack, E.W., E.A. Jacobson, S. Shipp and J.B. Anderson, Late 
Pleistocene-Holocene retreat of the West Antarctic Ice-Sheet 
system in the Ross Sea: Part 2 - Sedimentologic and 
Stratigraphic Signature, Geological Society of America 
Bulletin , 111, 1517-1536, 1999. 

Fastook, J.L., West Antarctica, the sea-level controlled marine 
instability: past and future, in Climate Processes and Climate 
Sensitivity, Geophys. Mono, vol 29, edited by J.E. Hansen and 
T Takahashi, pp 275-287, AGU, Washington, D.C., 1984 

Ferrigno, J.G., B.K. Lucchitta, K.F. Mullins, A.L. Allison, R.J. 
Allen, and W.G. Gould, Velocity measurements and changes 
in the position of Thwaites Glacier/iceberg tongue from aerial 
photography, Landsat images and NOAA AVHRR data, Ann 
Glaciol, 17,239-244, 1993. 

Giovinetto, M.B., The drainage systems of Antarctica: 
Accumulation, in Antarctic Snow and Ice Studies , Antarct 
Res. Ser. vol 2, edited by M. Mellor, pp 127-155, AGU, 
Washington, D C., 1964. 

Giovinetto, M.B. and C.R. Bentley, Surface balance in ice 
drainage systems of Antarctica, Antarct J. US., 20, 6-13, 
1985. 



VAUGHAN ET AL.: PINE ISLAND GLACIER 


255 


Goldstein, M., H. Engelhardt, B. Kamb, and R.M. Frolich, 
Satellite radar interferometry for monitoring ice sheet motion: 
application to an Antarctic ice steam, Science, 262, 1525- 
1530,1993. 

Hellmer, H.H., S.S. Jacobs, and A. Jenkins, Ocean erosion of a 
floating Antarctic Glacier in the Amundsen Sea, in Ocean, Ice 
and Atmosphere: Interactions at the Antarctic Continental 
Margin, Antarct. Res. Ser. vol 75, edited by S.S. Jacobs and 
R.F. Weiss, pp 83-100, AGU, Washington, D C., 1998. 

Hindmarsh, R.C.A., Qualitative dynamics of marine ice sheets, in 
Ice in the Climate System, NATO ASI. Series vol I 12, edited 
by W.R. Peltier, pp 68-99, Springer-Verlag, Berlin 
Heidelberg, 1993. 

Hughes, T.J., The weak underbelly of the West Antarctic Ice 
Sheet,-/ Glaciol, 27, 518-525, 1980. 

Jacobs, S.S., H.H. Hellmer, and A. Jenkins, Antarctic ice sheet 
melting in the Southeast Pacific, Geophys. Res. Lett., 23, 957- 

960, 1996. 

Jacobs, S.S. and J.C. Comiso, Climate variability in the 
Amundsen and Bellingshausen Seas, J. Clim., 10, 697-709, 
1997 

Jankowski, E.J. and D.J. Drewry, The structure of West 
Antarctica from geophysical studies. Nature, 291, 17-21, 
1981. 

Jenkins, A., D.G. Vaughan, S.S. Jacobs, H.H. Hellmer, and J.R. 
Keys, Glaciological and oceanographic evidence of high melt 
rates beneath Pine island glacier, West, Antarctica, J. Glaciol. , 
43, 114-121, 1997. 

Jezek, K.C., H.G. Sohn and K.F. Noltimier, The Radarsat 
Antarctic Mapping Project, IGARSS '98 - 1998 International 
Geoscience and Remote Sensing Symposium, Proceedings, 
Vol 1-5, Chapter 888, 2462-2464, 1998. 

Joughin, I., L. Gray, R. Bindschadler, S. Price, D. Morse, C. 
Hulbe, K. Matter and C. Werner, Tributaries of West 
Antarctic Ice Streams Revealed by RADARSAT 
Interferometry, Science, 286, 283-286, 1999. 

Jonas, M. and D.G. Vaughan, ERS-1 SAR mosaic of Filchner- 
Ronne-Schelfeis, Filchner Ronne Ice Shelf Programme 
Reports, 10, edited by H. Oerter, 47-49, AWI, Bremerhaven, 
Germany, 1996. 

Jones D.A. and I. Simmonds, A climatology of Southern 
Hemisphere extratropicai cyclones, Climate Dynamics , 9, 
135-145, 1993. 

Kellogg, T.B., R.S. Truesdale, and L.E. Osterman, Late 
quaternary extent of the West Antarctic Ice Sheet: new 
evidence from Ross Sea cores, Geology, 7, 249-253, 1979. 

Kellogg, T.B., D.E. Kellogg, and T.J. Hughes, Amundsen Sea 
sediment coring, Antarct. J US , 20, 79-81, 1985. 

Kellogg, T.B. and D.E. Kellogg, Recent glacial history and rapid 
ice stream retreat in the Amundsen Sea., / Geophys. Res , 92, 
8859-8864, 1987. 

Kellogg T B, D.E. Kellogg, E D. Waddington, and J.S. Walder., 
Late Quaternary' deglaciation of the Amundsen Sea. 
implications for ice sheet modelling, in The physical basis of 
ice sheet modelling, IAHS Pub 1 70 , edited by E D. 
Waddington and J.S. Walder, pp 349-357,Int. Assoc, of 
Hydrol. Sci., Wallingford, England, 1987. 


Legresy, B. and F. Remy, Altimetric observations of surface 
characteristics of the Antarctic Ice Sheet, J. Glaciol. , 43, 265- 

275, 1997.. _ 

Licht, K.J., N.W, Dunbar, J.T. Andrews, and A.E. Jennings, 
Distinguishing subglacial till and glacial marine diamictons in 
the western Ross Sea, Antarctica: implications for a last 
glacial maximum grounding line, Bulletin of the Geological 
Society of America, 111,91-103, 1 999 
Lindstrom, D. and D. Tyler, Preliminary results of Pine Island 
and Thwaites Glaciers Study, Antarct. J. U.S., 19, 53-55, 

1984. . „ . . 

Liu, H., K. Jezek, and B. Li, Development of an Antarctic digital 
elevation model by integrating cartographic and remotely 
sensed data: A geographic information system based 
approach.,/. Geophys. Res., 104,99-23,213, 1999. 

Lucchitta, B., C. Rosanova, and K. Mullins, Velocities of Pine 
Island and Thwaites Glacier, West Antarctica, from ERS-l 
SAR images, Ann. Glaciol., 21, 277-283, 1995. 

Lucchitta B.K., C.E. Smith, J. Bowell, and K.F. Mullins, 
Velocities and mass balance of Pine Island Glacier, West 
Antarctica, derived from ERS1-SAR. ESA Publication SP- 
361, 147-151, 1994. 

Lucchitta B.K., and C.E. Rosanova, Velocities of Pine Island and 
Thwaites Glaciers, West Antarctica, from ERS 1 -SAR images. 
ESA Publication SP-414 , 819-824, 1997. 

MacAyeal, D.R., Irregular oscillations on the West Antarctic Ice 
Sheet, Nature, 359, 29-32, 1992. 

Mantripp, D.R., J. Sievers, H. Bennat, C.S.M. Doake, K. 
Heidland, J. Idhe, M. Jonas, B. Reidel, A V. Robinson, R 
Scharroo, H.W. Shenke, U. Shirmer, F. Stefani, D.G 
Vaughan and D.J. Wingham, Topographic map (satellite 
image map), Filchner-Ronne-Shelfeis (2nd Edition) Map at 
1 2 000 000 , Institut fllr Angewandte Geodasie, Frankfurt 

am Main, Germany, 1996. 

Mercer, J.H., West Antarctic Ice Sheet and C0 2 greenhouse 
effect: A threat of disaster, Nature, 271, 321-325, 1978. 
Molnia, B.F. and P.R. Carlson, Surface sedimentary units of 
northern Gulf of Alaska continental shelf, Bulletin of the 
American Association of Petroleum Geologists, 62, 633-643, 
1999. 

Paterson, W.S.B., The physics of glaciers, Elsevier Science, 
Oxford, 480pp, 1994. 

Payne, A.J., A thermomechanical model of ice flow in West 
Antarctica. Climate Dynamics , 15, 115-125, 1999. 

Rignot, E.J ., Fast recession of a West Antarctic Glacier, Science , 
281,549-551,1998. 

SCAR, Antarctic digital database user's guide and reference 
manual , Scientific Committee on Antarctic Research, 
Cambridge, xi+156pp, 1993. 

Shimizu, H., Glaciological studies in West Antarctica, in 
Antarctic Snow and Ice Studies, Antarct Res Ser. vol 2, 
edited by M. Mellor, pp 37-64, AGU, Washington, D C., 
1964. 

Shipp, S., J. Anderson and E. Domack, Late Pleistocene- 
Holocene retreat of the West Antarctic Ice-Sheet system in the 
Ross Sea: Part 1 - Geophysical Results, Geological Society of 
America Bulletin , 111, 1486-1516, 1999. 



256 


THE WEST ANTARCTIC ICE SHEET: BEHAVIOR AND ENVIRONMENT 


Smith, A.M., The use of tiltmeters to study the dynamics of 
Antarctic ice-shelf grounding lines, J Glaciol , 37, 51-58, 
1991. 

Stenoien M.D., Interferometric SAR observations of the Pine 
Island Glacier catchment area. Unpublished Ph D. thesis. 
University of Wisconsin-Mad ison, 127pp, 1998. 

Stuiver, M., G.H. Denton, T.J. Hughes and J.L. Fastook, History 
of the marine ice sheet in West Antarctica during the last 
glaciation: a working hypothesis, in The Last Great Ice 
Sheets, edited by G.H. Denton and T.J. Hughes, Wiley - 
Interscience, New York, 1981. 

Stephenson, S.N., Glacier flexure and the position of grounding 
lines: Measurements by tiltmeter on Rutford Ice Stream, 
Antarctica, Ann. Glaciol . , 5, 165-169, 1984. 

Swithinbank C.W.M., Ice streams. Polar Record , 7, 185-186, 
1954. 

Thomas, R.H., Ice sheet margins and ice shelves. In: J.E. Hansen 
and T. Takahashi (Eds), Climate Processes and Climate 
Sensitivity, in Climate Processes and. Climate Sensitivity, 
Geophys. Mono, vol 29, edited by J.E. Hansen and T. 
Takahashi, pp 265-274, AGU, Washington, D C., 1984. 

USGS, Antarctica Sketch Map, Thurston Island - Jones 
Mountains , 1:500,000, USGS, Washington, D.C., 1993. 

Vaughan, D.G., Investigating tidal flexure on an ice shelf using 
kinematic GPS, Ann Glaciol , 20, 372-376, 1994. 

Vaughan, D.G., J.L. Bamber, M. Giovinetto, J. Russell, and 
A P R. Cooper, Reassessment of net surface mass balance in 
Antarctica, J. Clim ., 12, 933-946, 1999. 

Vaughan D.G. and J.L. Bamber, Drainage basin analysis and 
improved calculation of balance fluxes for West Antarctic ice 
streams and glaciers, in Abstracts, AGU Chapman 
Conference on the West Antarctic Ice Sheet, Orono, 
Sept. 1998, 1998. 


Weertman, J., Stability of the junction of an ice sheet and an ice 
shelf, J. Glaciol., 13,3-11, 1974. 

Whillans, I.M. and C. J. Merry, Ice-flow features on Ice Stream B, 
Antarctica, revealed by SPOT HRV imagery, J. Glaciol . , 39, 
515-527, 1993. 

Whillans, I.M., M. Jackson and Y-H Tseng, Velocity pattern in 
a transect across Ice Stream B, Antarctica, J Glaciol , 39, 
562-572, 1993. 

Williams, R.S., J.G. Ferrigno, T.M. Kent, and J.W. Schoonmaker, 
Landsat images and mosaics of Antarctica for mapping and 
glaciological studies, Ann. Glaciol, 3, 321-326, 1982. 

Wingham, D.J., A.J. Ridout, R. Scharroo, R.J. Arthem, and C.K. 
Schum, Antarctic elevation change from 1992 to 1996, 
Science, 282,456-458, 1998. 


C. R. Bentley and M.D. Stenoien, Geophysical and Polar 
Research Center, University of Wisconsin - Madison, Madison, 
Wisconsin 53706, USA. (email: bentley@geoIogy.wisc.edu) 

S. S. Jacobs, Lamont-Doherty Earth Observatory, University 
of Columbia, Palisades, New York 10964, USA. (email: 
sjacobs@ldeo.columbia.edu) 

T. B. Kellogg, Institute for Quaternary Studies, University of 
Maine, Orono, Maine 04469, USA. (email: 
tomk@iceage.umeqs.maine.edu) 

B.K. Lucchitta, U.S. Geological Survey, Flagstaff, Arizona 
86001, USA. (email: blucchitta@flagmail.wr.usgs.gov) 

E. Rignot, Jet Propulsion Laboratory, Pasadena, California 
91109, USA. (email: eric@adeire.jpl.nasa.gov) 

D. G. Vaughan, A.M. Smith, H, F. J. Corr, A. Jenkins, British 
Antarctic Survey, Natural Environment Research Council, 
Cambridge CB3 OET, U K. (email: d vaughan@bas.ac.uk) 



